Seismo Ionospheric Anomalies around and over the Epicenters of Pakistan Earthquakes

Abstract

Global Navigation Satellite System (GNSS)-based ionospheric anomalies are nowadays used to identify a possible earthquake (EQ) precursor and hence a new research topic in seismic studies. The current study also aims to provide an investigation of ionospheric anomalies associated to EQs. In order to study possible pre-and post-seismic perturbations during the preparation phase of large-magnitude EQs, statistical and machine learning algorithms are applied to Total Electron Content (TEC) from the Global Positioning System (GPS) and Global Ionosphere Maps (GIMs). We observed TEC perturbation from the Sukkur (27.8° N, 68.9° E) GNSS station near the epicenter of Mw 5.4 Mirpur EQ within 5–10 days before the main shock day by implementing machine learning and statistical analysis. However, no TEC anomaly occurred in GIM-TEC over the Mirpur EQ epicenter. Furthermore, machine learning and statistical techniques are also implemented on GIM TEC data before and after the M_w 7.7 Awaran, where TEC anomalies can be clearly seen within 5–10 days before the seismic day and the subsequent rise in TEC during the 2 days after the main shock. These variations are also evident in GIM maps over the Awaran EQ epicenter. The findings point towards a large emission of EQ energy before and after the main shock during quiet storm days, which aid in the development of lithosphere ionosphere coupling. However, the entire analysis can be expanded to more satellite and ground-based measurements in Pakistan and other countries to reveal the pattern of air ionization from the epicenter through the atmosphere to the ionosphere.

1. Introduction

Seismologists have found signs of EQ precursors in the form of changes in animal behavior, the groundwater level and temperature, electromagnetic signals, chemical emissions, and precursory fault slides and temporal differences in velocity, seismic activation, or precursory quiescence [1]. Furthermore, numerous seismological, geodetic, and other geophysical precursors, such as geoelectric, geomagnetic, geodetic, gravity precursors, hydro-seismology, or ground fluid, can be observed and monitored to improve EQ prediction procedures [2,3]. The Earth’s crust is complex and EQs are chaotic; therefore, EQ forecasts can be made statistically, and these statistical measures can be used to anticipate future EQs [4]. Moreover, the nature of the environment preceding the occurrence of EQs is challenging (e.g., EQ on land or in the ocean). On the other hand, the EQ precursors can be geophysical, geological, or atmospheric, even though EQs and the mechanisms that are at the root of these variations have yet to be fully defined [5,6].

As the peroxy bonds break, electrically charged carriers are released in order to create positive holes (p-holes) [7]. These p-holes can easily transport within the free volumes from trapped spaces inside the Earth. In a battery-like system, the electric current generated by each positive current transition results in low-frequency electromagnetic radiation and magnetic field changes.

Ionized molecules form amid the ground and air interface after p-holes arrive at the Earth. Thus, these molecules further rise to ionize the air as a probable cause of ionosphere anomalies [7,8]. Additionally, the possible relationship between the EQ activity and the changes in electron density is analyzed in different models [9,10]. The first theory is based on the piezoelectric effect that narrates the electron emission out of the trapped rocks. On the air-contact side of the granite block, ionization of molecules was observed after pressure was applied to its upper surface in a laboratory setting [11]. Furthermore, they showed that electrons and p-holes pass through the rock under pressure to the non-pressured side of the granite, where current flows between the pressure-applied side and the ending side of the rock block. Through the battery, this difference is reflected in a current flow. The negative pole of the battery cell is the pushed side, which is also loaded with positively charged molecules. The charge carrier keeps p-holes on the surface. The verification and confining variability in the Rocks–Earth surface are created by expanding spatially the environment of the experiment in the lithosphere–atmosphere–ionosphere system; however, prior to the commencement of a strong EQ, the ionization occurs in rocks where a fault ruptures due to stress. The system also distributes ionized air to the ionosphere layer, which in turn produces a higher electron density [12,13,14]. With this approach, EQs are classified by their magnitude with TEC anomalies. Ionospheric anomalies occur before EQs in the form of both positive and negative perturbations. P-holes are thought to cause positive ionospheric anomalies when the drift towards the ionosphere from an EQ induced an electric field and vice versa [14,15,16]. During adiabatically fluctuating conditions, positive ions are combined near the atmosphere and then move upward to the ionosphere. In the end, these complex reactions result in deviations in electromagnetic signals traveling through the ionosphere, which are detectable as ionospheric anomalies.

Electrons in the F-layer travel into lower layers because of the polarity of the pre-seismic electric field. This results in anomalous ionospheric parameters. Besides transporting electric fields to the upper atmosphere and eventually to the ionosphere, a thin layer of particles forms before EQs as a result of the ground’s emitted p-holes. A perpendicular electric field is created by transforming the vertical electric field on the ground surface into one that is perpendicular to the geomagnetic field lines. As a result of the zonal component, irregularities in the plasma density are observed before the EQ [10,17].

Research on ionospheric phenomena and precise positioning has been completed using dual-frequency global navigation satellite system (GNSS) receiver data. The total electron content (TEC) can be determined with dual-frequency receivers, providing users with additional information about the upper atmosphere. Ground and satellite measurements of TEC, foF2, and other ionospheric and atmospheric parameters were used to explore pre-and post-seismic anomalies triggered by large-magnitude EQs [18,19,20,21,22,23,24]. Through various models, more theoretical reports were undertaken to elucidate the fundamental mechanism of seismo-ionospheric anomalies (SIAs) propagation. In addition to SIAs, the effects of a geomagnetic storm must be considered while examining ionosphere disturbances [15,17]. On account of the seismic preparation period, several researchers discussed various strategies for detecting and monitoring the co-located and synchronized abnormalities from both ground and space observations above the epicenter of an upcoming major shock. For example, for the Chi-Chi event in Taiwan in 1999, a significant link was established statistically between the EQ and ionospheric abnormalities [25]. Shah and Jin [24] used TEC to identify ionospheric abnormalities interrelated with several terrestrial EQs stirring at various fault systems, accentuating the significance of monitoring SIAs associated with Mw > 6.0. Oyama et al. [26] found an unusual fall in ion density data from the US satellite DE-2, which they linked to Mw 7.5. Ionosphere differences at different altitudes are most likely explained by the creation of electric fields around the upcoming EQ’s epicenter.

The ionosphere is a complex layer in the upper atmosphere that fluctuates in space and time, owing to geomagnetic activity and ionization by solar extreme ultraviolet (EUV) radiation. Solar and geomagnetic factors are typically assigned to all energy input for ionosphere fluctuations, but other influences have been largely overlooked [10,27]. Recent studies have shown the importance of ionospheric disturbances for EQ precursors monitoring [28], and major studies have focused on ionospheric perturbations due to EQs [10,29,30,31]. Additionally, ionospheric anomalies are generated by several mechanisms, and three of them have received much attention, among them, the Freund [32], proposed p-hole activation due to stress build up, and Ondoh [33], theorized radon emanation-based model, while the Namgaladze et al. [17] model of the atmospheric gravity wave received attention from all over the world.

A variety of methods are available for studying seismo-ionospheric abnormalities. The noteworthy approaches are statistical-based seismo-ionospheric analysis [30,34,35], case study analysis [26], and physical model-based analyses for the elucidation of ionospheric anomalies [12]. Statistical methods can be used to show a link between the EQ and the anomalies it causes. On the other hand, the use of case studies and physical model analysis is useful when explaining the physical framework of contemporary seismo-ionospheric anomalies. In accordance with the previous discussion, seismic precursors are still a subject of debate, and there was no dedicated satellite mission to monitor EQ-induced ionospheric irregularities prior to the DEMETER satellite [36].

Traditional anomaly detection approaches are based on statistical analysis and rely heavily on past data. Any statistical methods are incapable of replicating past observations, and patterns are far from the forecasting of EQ anomalies. Artificial Neural Networks (ANNs) have been widely used as time series forecasters to show the difference between the predicted and observed values at a predefined threshold value in regard to EQ anomalies [37]. In this study, we investigate the ionospheric anomalies associated with two Pakistani EQs, using statistical and machine learning approaches for the seismic events of Mw > 5.3 in the northern and southern parts of the country.

2. Materials and Methods

2.1. Materials

In this study, we analyzed two seismic events of Mw > 5.3 in the northern and southern parts of Pakistan. The Awaran EQ of Mw 7.7, which has a 15 km hypocentral depth, occurred on 24 September 2013 (Figure 1). Moreover, the epicenter was located at 113 km northwest of Bela and 61 km northeast of Awaran on the USGS website http://EQ.usgs.gov/EQs/eventpage/usb000jyiv#general_summary (accessed on 08 January 2023). The EQ epicenter has Latitude: 26°09′60.00″ N and Longitude: 65°29′59.99″ E. The maximum felt intensity was VIII (severe) on the Modified Mercalli scale, as is on the record of the United States Geological Survey (USGS). About a minute after the main shock, panic in the region resulted in the loss of 825 lives and significant infrastructure damage. The EQ was followed by a series of strong aftershocks, with a major of M_w 6.8 on 28 September 2013. The affected areas were Awaran, Tirtej, Gashkore, Nok Jo, Parwar, Dandar, and Hoshab. The EQ occurred due to an oblique-strike slip type motion. The tremor occurred in all Pakistan and its neighboring countries: Iran, Afghanistan, and India. A cumulative moment of 5.4 × 10²⁰ Nm was released by this rupture.

Another case study includes the Mirpur EQ of shallow hypocentral depth of 10 km, which hit the Pakistani region on 24 September 2019 (Figure 1). Its epicenter was located at 7 km SSE of Mirpur, Pakistan (33.1° N, 73.7° E). The maximum felt intensity of VII (very strong) on the Modified Mercalli scale was reported by USGS. There was serious damage to the lives and the infrastructure of the people in Mirpur District, resulting in 40 deaths and more than 850 injured.

Furthermore, a number of bridges and many roads were demolished, especially within 14 km of the main seismic zone. The epicenter of the EQ was found 3 km south of Mirpur and 20 km north of Jhelum in Pakistan. More information about the EQ is available on the USGS website (https://EQ.usgs.gov/EQs/eventpage/us60005mqp/executive) (accessed on 8 January 2023). We obtained the GPS TEC data from the Sukkur GNSS station in Pakistan for the Mirpur EQ. Similarly, the ionospheric anomalies are also studied from the GIMs of IGS for the Awaran, Pakistan EQ in 2013. In this study, the Slant TEC (STEC) from the Sukkur GNSS station is estimated as the number of free electrons in a square meter along the line of sight between the satellite and the receiver, as given by the equations below.

$$STEC=\frac{f_{1 }^2{f}_2^2 }{40.28\left(f_1^2-f_2^2\right)}\left(L_1-L_2+{\lambda }_1\left(N_1+b_1\right)-{\lambda }_2\left(N_2+b_2\right)+ϵ\right)$$

(1)

$$STEC=\frac{f_{1 }^2{f}_2^2 }{40.28\left(f_1^2-f_2^2\right)}\left(P_1-P_2-\left(d_1-d_2\right)+ϵ\right)$$

(2)

In the above equations, the carrier phase frequencies are presented by f₁ and f₂, the pseudo-range is denoted as L, the delay path of the signal of the carrier phase observations is P, the signal wavelength is λ, and the ray path uncertainty is N. Here, d and b denote the biases of the consequent signal pseudo-range and instrumental carrier phase, and ϵ is the random error in the signal. The STEC is converted to VTEC using the following equation [8].

$$VTEC=\mathrm{STEC} \times \cos \left(arcsine\left(\frac{RsinZ}{R+H}\right)\right)$$

(3)

In this equation, Z is the elevation angle of the satellite, and R and H are the Earth’s radius and the ionosphere height, respectively.

2.2. Methodology

The anomalous VTEC values for the EQs in the northern and southern parts of Pakistan are further investigated in a statistical analysis. In order to detect abrupt changes in VTEC derived from GPS satellites, confidence bounds are applied using the following equation.

$$Confidence Bounds=u\left(TEC\right)\pm 1.5 \sigma \left(TEC\right)$$

(4)

The $u\left(TEC\right)$ and $\sigma \left(TEC\right)$ are mean and standard deviation of the observed values, respectively. In order to detect anomalies, we observed abnormal VTEC values out of the bounds. The deviation from the normal distribution of VTEC is indicated by a value exceeding the upper bound or below from the lower bound.

$$Anomalous Behavior=\left\{\begin{array}{c}\mathrm{Anomaly} TEC>UB\\ \mathrm{No} \mathrm{Anomaly} UB<TEC<LB\\ \mathrm{Anomaly} TEC<LB\end{array}\right.$$

(5)

A class of intelligent systems called Artificial Neural Networks can discover patterns based on a few priori assumptions and examine complicated functional relationships in data to predict a phenomenon. Despite unknown or complex underlying laws governing a time series, neural networks can capture the autocorrelation structure in the series [38]. In a large number of applications, neural networks have been very successful in modelling and forecasting nonlinear time series [39]. ANNs use simple processing units, called neurons, to store data related to the learning process. They are connected in a network by numerous connections, and each stores a certain amount of data. In ANN, the product x_j·w_ij is produced by multiplying a strength x_j by a weight w_ij. A neural network’s output is represented by the following: y_i = f(x_jw_ij), where i and j are indexes of a neuron in the hidden layer and input, respectively [38].

Multi-Layer perceptron is a feed-forward network, in which only a forward connection exists, and neurons are grouped together. It consists of n input units, along with k hidden units, but has only one output unit. Connection weights are described as Wij and Tj, which are the weight from the ith input unit to the jth hidden unit and the jth hidden unit to the output unit, respectively. The most important part is the hidden layer having nodes in connection with the input, as well as the output, layer. For anomaly detection, complete data are split into two sets, training and test, where the training set is used for the construction of the neural network, while the test set is utilized for the assessment of the predictive error, whereas the network’s connection weight is determined by the training process. Measured values may be considered anomalous if the prediction error exceeds the predefined threshold. To assess the minimum error while training, the learning algorithm, so as the activation function, was obtained by the process of iteration. When configuring the best suitable network, premium significance was given to the parameters that directly impact the predictive error, encompassing the pattern input number, hidden layer, lag value, etc. With regards to the output layer activation function, the linear function gave the best results. A tan sigmoid function is used to transfer the information for all hidden nodes. Sigmoid functions can be thought of as a kind of mathematical function:

$$f\left(x\right)=\frac{2}{1+e^{-2z}}-1$$

(6)

Training was carried out using the well-known Levenberg–Marquardt algorithm. Predictions were made using N observations x₁, x₂, …, x_N as the training data, and the remaining observations x_N+1, x_N+2, …, x_N+m as the test data. The time series underlying pattern is discovered using lagged observations that correspond to the input nodes (p). The network capability of the prediction or learning is highly affected by the various input nodes [38]. For our analysis, the number of nodes for the input, hidden, and output layers was three, two, and one, respectively. Specifically, each training set contains four observations that constitute one pattern vector, of which three are inputs and one is the output value.

The proposed network has the following training patterns:

X_N = f(X_N−3, X_N−2, X_N−1); N = 1, 2, 3….

(7)

During the training process, the connection weights are optimized to minimize the prediction error (PE). Here is an equation for PE:

$$PE={{\displaystyle \sum }}_{i=4}^N{X}_i-{\hat{X}}_i$$

(8)

${\hat{X}}_i$ is the network’s output.

The patterns for testing are:

X_N+m = f (X_N+m−3, X_N+m−2, X_N+m−1)

(9)

We use recurrent neural networks to evaluate anomalies in the VTEC time series, using nonlinear autoregressive models with exogenous inputs (NARX) [40]. This NARX model for the time series leverages regressing to estimate the dependent output signal y using previous estimations of both the output and independent input signals (t).

The equation of the NARX model is showed below.

$$\mathrm{y}\left(\mathrm{t}\right)=\mathrm{f}(\mathrm{y} (\mathrm{t}-1), \mathrm{y}(\mathrm{t}-2), \dots , \mathrm{y}(\mathrm{t}-{\mathrm{n}}_{\mathrm{y}}), \mathrm{x}(\mathrm{t}-1), \mathrm{x}(\mathrm{t}-2), \dots , \mathrm{x}(\mathrm{t}-{\mathrm{n}}_{\mathrm{x}}\left)\right)$$

(10)

In the training algorithm, Levenberg–Marquardt optimization is used. Additionally, during implementation of the prediction procedure, N observations consisting of y₁, y₂, …, y_N and y_N+1, y_N+2, …, y_N+m are selected as training and test sets, respectively. We can identify the underlying pattern of a time series by counting the lagged observations corresponding to input nodes. Neurons receive observed values from a variety of sources and feed them into the input layer. In this study, the corrected LST values, their accompanying times, and the deviation from the mean LST values are inputs. All neurons in the hidden layer are affected by the weights collected from the previous input layer. An effective combination of hidden neurons engages the transfer function.

The network’s training patterns are as mentioned below.

y_N = f(y_N−3, y_N−2, y_N−1, t_N−3, t_N−2, t_N−1)

(11)

Prediction performance assessment involves finding optimal weights to minimize the prediction error (PE).

The PE equation can be inscribed as:

$$\mathrm{PE}={{\displaystyle \sum }}_{\mathrm{k}-0}^{\mathrm{N}} \left(\widehat{y}\left(\mathrm{t}-\mathrm{k}\right)-y\left(\mathrm{t}-\mathrm{k}\right)\right)$$

(12)

$\widehat{y}$ is the network’s output.

The testing patterns are:

$${\mathrm{y}}_{\mathrm{N}+\mathrm{m}}=\mathrm{f}\left({\mathrm{y}}_{\mathrm{N}+\mathrm{m}-3},{\mathrm{y}}_{\mathrm{N}+\mathrm{m}-1},{\mathrm{y}}_{\mathrm{N}+\mathrm{m}-1},{\mathrm{t}}_{\mathrm{N}+\mathrm{m}-3},{\mathrm{t}}_{\mathrm{N}+\mathrm{m}-2},{\mathrm{t}}_{\mathrm{N}+\mathrm{m}-1}\right)$$

(13)

Finally, the anomaly is identified by determining whether the difference between the expected and real data exceeds predefined limits.

3. Results

In this study, TEC data above the epicenter for possible EQ anomalies are investigated using confidence bounds, NARX, and MLP techniques within the period before and after the main shock day. The detail results for both the case studies is mentioned below.

3.1. Case Study I

The significant ionospheric TEC anomaly was observed on 22 September 2019, specifically during LT = 10–12 h beyond the bounds (Figure 2), while no irregular value was observed after the main shock in the form of a post-EQ anomaly. Hence, we can say that only pre-seismic ionospheric anomalies were noticed in the case of the Mirpur EQ. The anomaly on 22 September 2019 is also investigated in the ionospheric maps over the epicenter in IGS-retrieved GIMs (Figure 3). One can see no prominent TEC cloud in GIMs around the epicenter, which may be due to the lower magnitude of this EQ. On the other hand, it is also possible that the seismic energy may not have reached ionospheric heights. The GIM maps also have low resolution bi-hourly over the epicenter, and it is also possible that the low-resolution ionospheric GIM maps may also have not been capable of detecting short-term ionospheric perturbations. A series of short-term ionospheric anomalies of about <1 TECU above the upper confidence bound occurred in the NARX analysis within 6 days before till 3 days after the main event (Figure 4).

However, the NARX method showed a small anomaly on 20 September 2019 for the few final hours of the day around 3.6 TECU, and on 22 September, a small anomaly of 2 TECU was detected. Since the original VTEC from the GPS only has an anomaly on 22 September, the other small peaks in NARX may be the outliers. NARX was used to enhance the peaks in time series forecasting. Similarly, the MLP technique observed an obvious deviation on 22 September 2019 of the same magnitude as the original GPS VTEC. The deviated value beyond the upper bound from MLP is 7.14 TECU (Figure 5), whereas the GPS value is 10.65 TECU in Figure 2. All the anomalous days before and after the main shock of the Mirpur EQ are showed in

Table 1. Brief description of seismo-ionosphere anomalies.

Method	Pre-EQ Anomalies				Post-EQ Anomalies
	Mirpur		Awaran		Mirpur		Awaran
	Day	Deviation	Day	Deviation	Day	Deviation	Day	Deviation
IQR	−2	2 TECU	EQ day	4 TECU	Nil	Nil	Nil	Nil
NARX	−5, −2	4 TECU	EQ day	5 TECU	1–2	2 TECU	Nil	Nil
MLP	−2	3.5 TECU	EQ day	7 TECU	Nil	Nil	1 & 3	3 TECU

.

3.2. CASE Study II

One can see clear ionospheric anomalies in the GPS VTEC of the Sukkur GNSS station on the day of the main shock, which justifies the seismic energy emanation to the ionosphere (Figure 6). The anomalous pattern is observed in the analysis of VTEC from GPS for 10 days before and 5 days after the main shock day, which is studied for the Awaran EQ. While analyzing the temporal data of TEC, we observed a 0.67 TECU rise in the electron cloud beyond the upper confidence bound on the EQ day. This variation in TEC is further studied in reconstructed GIMs over the epicenter of the Awaran EQ on 24 September 2013 (Figure 7). It is interesting to see anomalous TEC clouds of > 5 TECU over the epicenter during LT~10–12 h on the EQ day. The anomalous TEC clouds initiated during LT~9 h circulated over the epicenter till LT~12 h and vanished after LT~12 h (Figure 7). These ionospheric clouds over the epicenter of the Awaran EQ endorsed the temporal VTEC variations from the Sukkur GNSS station. It is also different from the Mirpur EQ, whose temporal and spatial TEC showed anomalous TEC concentrations on the EQ day. These kinds of ionospheric anomalies several hours to several days after the main shock were widely reported [37], but interestingly, we also logged a few post-seismic irregularities in TEC. The detail about the ionospheric anomaly before/after the main shock can be found in Shah et al. [21]. The irregularities in TEC before the main shock support the fact that the propagation of EQ energy via the atmosphere to the ionosphere is validated under the assumption of lithosphere–atmosphere–ionosphere coupling [18,33,41].

Nevertheless, the observed irregularity on the day of the EQ, i.e., 24 September 2013, was 5 TECU during the afternoon hours, i.e., LT~11–12 h in the NARX-predicted VTEC (Figure 8). While evaluating TEC data using Neural Networks, we observed the same EQ day ionospheric anomaly with a more significant peak. It is calculated as the difference (DX_i) between the actual value (X_i) and the predicted value $(\widehat{X_i})$ outside the defined bounds.

The implemented NARX model on the original VTEC within the upper and lower bounds endorsed the ionospheric anomaly on the EQ day for the Awaran event. The anomaly beyond the upper bound for MLP is 3.5 TECU on the main shock day as a seismic response (Figure 9). There also occurred a slight anomaly on 18 September 2013, and the deviation observed was 1.43 TECU. However, it did not cross the bounds to make it a SIA. Moreover, the machine learning techniques only show the ionospheric anomaly to a minute level in both cases and do not assist in finding a new VTEC precursor. All the anomalous days before and after the main shock of both the EQs, with their deviations from the upper/lower bound, are showed in Table 1.

The geomagnetic storm and solar activity indices showed no clear variation associated with both main shocks (Figure 10). The geomagnetic indices, including Kp, Dst, F10.7 and AE, predicted a quiet pattern of solar activity; thus, the ionospheric perturbations are legitimately induced by the EQ main shock event. Figure 10 depicts that Kp < 4 and Dst > −30 nT before the main shock day and also around the anomalous day, a similar pattern was observed, hence proving that the ionospheric anomalies were due to the EQ for both events.

4. Discussion

In this paper, we found significant variations in GPS TEC before the main shocks of the Awaran and Mirpur events within a 5–10 day window of the EQ preparation period (Figure 2 and Figure 6). These abnormal variations propagated to the ionosphere via a complex process through the atmosphere and lower ionosphere during the main shock preparation period before the EQs. This process of EQ and ionosphere coupling is clearly showed by Daneshvar and Freund [42] in the form of atmospheric blocking and seismic events in the Middle East region, with a correlation of 2.88 for establishing a coupling hypothesis between the lithosphere and atmosphere, followed by a clear variation in the ionosphere. The reason of the ionosphere variations prior to the main shocks is the drifting of abnormal atmospheric blocking toward the ionosphere, which caused delays in the GPS signals. The lithospheric variations prior to the EQ in the form of abnormal electric field generation, electric flux exchange, radon gas emission, and the p-hole hypothesis can clearly describe a background for the connection between the atmospheric and ionospheric precursors with seismic events, which is also clear in our case studies (Figure 3 and Figure 7). Previous reports highlighted the future main shock events as the source of coupling for the mechanisms of lithospheric activities, atmospheric variation, and abnormal ionospheric irregularities [43,44]. Therefore, we also hypothesized the development of ionosphere anomalies associated with the Pakistani events due to the p-holes’ emission. However, the ionospheric anomalies triggered due to radon emanation in the LAIC model cannot be omitted for the physical processes controlling the possible EQ and ionospheric precursory phenomena [13].

The abnormal ionospheric anomalies in the GPS data can be discussed under the hypothesis of two prominent physical processes before the main shocks during the seismic preparation period; one addressed the release of radon and other gases from the lithosphere before the EQ for producing air ionization (e.g., [45,46,47]), and the second involved p-holes’ emission from the stressed rocks as charge carriers at hypocentral volume to generate an electric field at the air–ground interface [30,39]. Before the EQ during the main shock preparation period, increasing strains triggered electronic charge carriers due to the build-up of tectonic stresses deep inside the Earth, which were active charge carriers (p-holes) within a week before the EQ day. The p-holes deficient in electrons are highly mobile and have significant spread ability due to the wave function of O^- in the matrix of O^2- within the rock columns of the sublattice of oxygen. As soon as these p-holes reach the Earth’s surface (lithosphere), they accumulate over the lithosphere at the ground–air interface at different topographic highs to produce electric fields. Furthermore, these electric fields over the epicentral regions are microscopic in nature, but they have the ability to reach and spread around millions of volts per centimeter. This abrupt electric field causes air ionization, mostly O² and many other molecules over the lithosphere in the nearby air. Over the EQ preparation region, the model of global electric current propagation, Earth electric charge carriers, and ionization of air raised to high altitude by the explanation of Freund et al. [48]. The high altitude of the seismic regions further intensifies the electric field due to its intrinsic nature [49], and the case studies in this paper also occurred in high and heterogeneous regions. Previous reports on various atmosphere and ionosphere precursors also highlighted the variations at different altitudes from satellite observation [50,51]. After the loss of an electron from O² to the ground state, it came with O²⁺ as positive ions in the atmosphere. The air with more and more O²⁺ in the seismic preparation period rose further upward as a result of the condensation of moist particles and latent heat release. The previously discussed moisture condensation further raises anomalous latent heat, and the air ionization causes an increase in air velocity and flow towards the upper atmosphere and ionosphere [52]. The enhancement in air buoyancy generates thermal up gradation, which aids in the transportation of ionized air to high altitudes over the EQ epicenters and culminates in the form of ionospheric anomalies.

During the EQ preparation period before the main shock, the hypothesis of the release of radon and other gases from the epicentral regions is also responsible for the air ionization, condensation of atmospheric moisture, and abnormal precipitation [1,42,43,44,45,46,47]. The reason is that radon increases the ionization rate in the atmosphere and nucleation [53]. This air-laden radon can further drift in the atmosphere towards the lower ionosphere, which is followed by an increase in the upper ionosphere. This is clear before the seismic events over the epicentral regions in the observations of different satellites [54,55,56,57]. The phenomena of radon decay explained the production of an equal number of positive and negative ions at the surface-to-air interface, followed by the recombination process after reaching a threshold value, thus causing a decrease in the air ionization [58,59,60,61,62]. Previous reports also analyzed these variations [63,64,65,66,67]. Furthermore, the number of condensation nuclei formed the air droplet release of latent heat and buoyancy. By contrast, the p-holes hypothesis can produce more positive airborne ions over the epicentral region in the air to form variations in the atmosphere. On the other hand, the electrostatic imbalance in the lower edge of the troposphere can also clear the variations of ionospheric irregularities due to seismic events. Furthermore, the satellite-based atmospheric and ionospheric anomalies provided deep insight on future main shock EQs [68,69,70,71,72]. These variations become clearer with a precise method from the satellite data [73,74]. For this, we need more analysis with multiple satellite observations to study the atmospheric and ionospheric variations. Further, more adequate methods can also assist in the delineation of appropriate seismic precursors in the form of pre-EQ situations.

5. Conclusions

In this study, we have investigated the ionospheric anomalies associated with Pakistan EQs using statistical and machine learning approaches for Mw > 5.3 in the northern and southern parts of the country. The following conclusions are made from the present work.

• Ionospheric anomalies occurred before the Mirpur EQ as pre-SIA within 5 days before the main shock, and variations on the Awaran EQ day occurred as a seismic response on the main shock day. Moreover, the intensity of the VTEC anomaly for the Awaran EQ was higher than the Mirpur event due to the magnitude difference.
• The differential GIM maps showed no clear electron cloud over the epicenter of the Mirpur EQ due to the low magnitude of the event. On the other hand, the Awaran EQ induced significant TEC clouds over the epicenter during LT~10–12 h.
• The most apparent anomalies occurred before the Mirpur EQ and on the EQ day of the Awaran event as a pre-SIA and seismic response, respectively. Similarly, no clear post-EQ anomalies occurred in the case of both EQs. These results suggest that abrupt seismic variations triggered by the EQs appeared in the form of the emanation of energy from the EQ-prone region to the ionosphere in the seismic preparation period. In conformity with previous studies and our conclusions, we believe that the seismo-ionospheric and thermal anomalies can positively contribute to the prediction of EQs.
• Machine learning can only assist in enlarging the peak and variations of the existing abnormal VTEC values and can no longer help in finding new VTEC precursors.
• The geomagnetic indices also observed no effect of storms for the same observation period before and after the EQs. Thus, we can conclude that the observed TEC anomalies were triggered by the seismic events.
• The GIMs provide a clear picture of authentic ionosphere anomalies in the case of the Mirpur EQ. They show that low-magnitude EQs cannot propagate energy to ionospheric heights. However, we can monitor the response of TEC to seismic events by integrating observations based on the worldwide impact of such natural phenomena, and it might be utilized as a useful tool to forecast possible seismic activity. Ionosphere-seismic studies are a continuing process and one of the primary drivers of ionosphere variability, and such consequences can be noticed at least a few days before an EQ.