*4.2. Time Series Analysis of Epicenter TEC*

To observe the affected hours more closely, TEC data were obtained for 15 days at the intersection of the four epicenters and two seismogenic zones, as shown in Figure 4, where the red curve is the upper quartile threshold, the green curve is the lower quartile threshold, the blue curve is the actual TEC value, and the pink and black vertical bars are the earthquake occurrence times. If the TEC data exceeded the upper and lower quartile thresholds, the TEC deviations were drawn on the red histogram at the bottom of the graph. TEC deviations can be seen beyond the upper limits, 3–4 days before successive earthquakes. We observed slight seismo-ionospheric anomalies 10 days before the earthquake. As shown in Figure 4a, at 07:00 UT on 16 November 2019 (10 days before the earthquake), TEC first appeared as a negative anomaly of −0.92 TECu, and at the same time of the next day (17 November), it also appeared and reached −0.93 TECu. Then, the negative anomaly gradually weakened and became a positive anomaly. The positive anomalies were initiated at 08:45 UT on 22 November (4 days before the earthquake), and these positive anomalies, occurred between 08:45 UT on 22 November and 15:45 UT on 23 November (3 days before the earthquake). The maximum time interval of the anomaly was only 1 h. The peak of the positive anomalies appeared at 17:45 UT on 21 November (5 days before the earthquake), and reached 2.54 TECu. The anomalous characteristics of the epicenters of the Bosnia and Herzegovina earthquakes were consistent with those of the Albania earthquake, as shown in Figure 4b. This phenomenon could be explained because the epicenters of the Bosnia and Herzegovina earthquakes were only 137.8 km apart from the epicenter of the Albania earthquake and the two seismogenic zones coincided. However, as shown in Figure 4c, with a single seismogenic zone at the epicenter, there was a negative anomaly of TEC = −2.91 TECu at 23:30 UT on 15 November (11 days before the earthquake), which is similar to the characteristics of the two previous epicenters, while the intensity of the seismo-ionospheric anomalies had no continuity. Then, positive anomalies appeared at 17:30 UT on 22 November (5 days before the earthquake), with the peak at 03:00 UT on 23 November (3 days before the earthquake), and reached 3.42 TECu. Figure 4d shows the seismic ionospheric anomalies in the area where the Albanian seismic zone meets the Greek seismic zone. The positive anomaly characteristics at the intersection of seismogenic zones were very similar to those in the Albania seismogenic zone, with both showing the characteristics of continuous anomalies. Compared with the single earthquake area, the seismic ionospheric anomalies were significantly different. In addition, earthquake-induced VTEC anomalies occurred more frequently within a 3–10-day window with the approach of the earthquake day. The intersection of seismogenic regions resulted in more frequent ionospheric fluctuations. These results suggest that there is a certain correlation between the intersection of seismogenic regions in time and space, and ionospheric disturbance.

#### *4.3. TEC Changes in the Balkan Peninsula Seismic Swarm*

In order to study whether CID are caused by various atmospheric waves (Rayleigh surface waves, and internal gravity waves) generated by the seismic swarm, we investigated the TEC responses to the Balkan Peninsula seismic swarm [28,29]. Figure 5 shows the VTEC time series and the trajectory of the ionospheric penetration point (IPP) observed by the GPS station in the earthquake area, where the orbit color matches the color of the observation station. The pentagram on the orbit is the location of the satellite when the earthquake occurs, and the pentagram color corresponds to the color of the earthquake center (when the red earthquake center occurs, it corresponds to the red five-star position on the IPP path, and the blue earthquake center corresponds to the blue five-star). The red arrow is the direction in which the satellite moves.

**Figure 5.** *Cont*.

**Figure 5.** The left panels present the time series of the original VTEC changes observed by different satellites at different stations and the right panels present the corresponding IPP trajectory of the satellite. Additionally, red, blue and green lines indicate the occurrence of an earthquake in the left figure (the color corresponds to the source in the right figure). The small stars on the track represent the IPP position of the satellite when the earthquake occurs (the color corresponds to the source in the figure) and the red thread indicates the moving direction of the satellite in the right figure. (**a**) The time series of the original VTEC changes observed by satellite 16 and the IPP trajectory of the satellite. (**b**) The time series of the original VTEC changes observed by satellite 27 and the IPP trajectory of the satellite. (**c**) The time series of the original VTEC changes observed by satellite 10 and the IPP trajectory of the satellite. (**d**) The time series of the original VTEC changes observed by satellite 08 and the IPP trajectory of the satellite. (**e**) The time series of the original VTEC changes observed by satellite 09 and the IPP trajectory of the satellite.

In the left plane of Figure 5a, we show the VTEC series from 1:00 UT to 5:00 UT recorded by satellite 16 visible from different stations (except for the ZADA station, the IPP longitude range is 11◦–13◦E) in the seismogenic area. For all stations, CID clearly appear after the earthquake, with time lags of 1–2 h, This phenomenon is consistent with the conclusion of Astafyeva [30,31], and the time needed for acoustic waves to travel from the surface to the IPP. Interestingly, we found that when we compared the VTEC sequence of ZADA station with that of TUC2 station, ZADA station is at the edge of the seismogenic area, IPP also does not pass over the seismogenic area, and the fluctuation amplitude of VTEC sequence is not very obvious. On the contrary, although the TUC2 station is not in the earthquake preparation area, the IPP track is located above the earthquake preparation area and close to the focal center, and it can be clearly seen that it has been strongly disturbed. In addition, it can also be observed that within 2 h before the arrival of the first main earthquake, all stations were disturbed to varying degrees. The perturbation frequency increased as the IPP positions got closer to the source center, with a fluctuation range of 0.1–0.5 TECu, and similar phenomena were observed on other satellites.

The TEC shows rising or falling curvature due to satellite elevation and latitude changes, so it is difficult to observe abnormal conditions. Nevertheless, we observed from the IPP path of satellite 27 (Figure 5b) and satellite 09 (Figure 5e) within 2 h, and its latitude change rate was relatively low. It is obvious that compared with satellite 27 (Figure 5b), when the 6.3 Mw Albanian earthquake struck, CID was observed at all stations on satellite 27. In particular, the ORID station is close to the center of the source, and the IPP is also closest to the center of the source, and its variation in amplitude is large. Figure 5e shows VTEC changes over the seismogenic area after three earthquakes. From the figure on the right, we can see that the IPP of observation satellite 09 (Figure 5e) and the POZE, ZADA and DUB2 stations monitor the VTEC in the upper half of the seismogenic area. In addition, satellite 09 was observed at other stations, and its IPP position was located in the lower half of the seismogenic area. The VTEC observed at the ORID, LARM, and PAT0 stations were disturbed to the same extent, and VTEC change curves showed similar changes. However, at the TUC2 station outside the seismogenic area, it was observed that the VTEC at the edge of the seismogenic area was significantly different from other stations at UT 10:00. Interestingly, the overall trend of VTEC still showed a downward trend, indicating that the decrease in TEC after the earthquake may result from electronic transport. Therefore, it can be inferred that the earthquake swarm will continuously interfere with the ionosphere over the seismogenic area, and this interference will continue to increase with the increase in the number of earthquakes. It is interesting that within 2 h after the earthquake, the VTEC over this area decreased sharply (in the upper half of the seismogenic area). In order to analyze this abnormal phenomenon, the earthquake swarm likely caused the formation of ionospheric holes in the region [32]. In Figure 6b, at UT 9:00–10:00, the VTEC (satellite 22) in the seismogenic area observed by the three stations was lower than 10 tecu. However, after UT 11:00, when satellite 09 moved into the region, VTEC still showed a downward trend, which preliminarily indicated that the ionosphere in this region was abnormal. Figure 6a shows the change in VTEC on the day before the earthquake. It can be seen that its change characteristics align with the daily change characteristics of the ionosphere. With time, VTEC gradually increases.

**Figure 6.** POZE, ZADA and TUB2 station observed the IPP position and VTEC changes of satellite 09 and satellite 22 at UT 9:00–12:30. The left panels present TEC disturbance (TEC) along the IPP trajectories (trajectories are from left to right) in the vicinity of the epicenter on the previous day on November 25 (**a**) and the right panels show the day of the event on November 26 (**b**).

#### *4.4. Observed Coseismic Ionospheric Disturbances in the Balkan Peninsula Seismic Swarm*

To further analyze the CID, we examined TEC data with elevation angles higher than 20◦. The TEC disturbance (dTEC) was obtained by employing high pass filtering.

In Figure 7, we compare the IPP trajectories and dTEC of satellite 16, 10, and 09 observed at the GPS stations. Abnormalities can be seen in the dTEC after the onset of the mainshock (26 November 2019) on earthquake day, but no such anomaly can be noted on the previous day or on the day after, respectively. In addition, it can be clearly seen that the anomaly of satellite 16 on the mainshock day is relatively weak, and the anomalies of satellite 10 and satellite 09 are more intense on the mainshock day. The IPPs of satellites 10 and 09 are approaching the source center. After three earthquakes in this area (25 November 2019), we noticed that the variation in the dtec amplitude of satellite 09 was about 0.15 tecu. Therefore, from this spatiotemporal analysis, we concluded that the intensified TEC disturbances in Figure 7 were likely of seismic origin.

(**a**)

**Figure 7.** *Cont*.

**Figure 7.** The upper panels present TEC disturbance (dTEC) was observed on three consecutive days (25–27 November 2019, event day on 26 November 2019). The lower panels show their corresponding IPP trajectory of the satellite. (**a**) TEC disturbance (dTEC) on the previous day on 25 November 2019, on the event day on 26 November 2019, and on the day after on 25 November 2019, and the IPP trajectories of the GPS Sat. 16. (**b**) TEC disturbance (dTEC) on the previous day on 25 November 2019, on the event day on 26 November 2019, and on the day after on 25 November 2019, and the IPP trajectories of the GPS Sat. 10. (**c**) TEC disturbance (dTEC) on the previous day on 25 November 2019, on the event day on 26 November 2019, and on the day after on 25 November 2019, and the IPP trajectories of the GPS Sat. 09.

#### **5. Discussions**

In this paper, the PEIAs associated with a seismic swarm in Balkan Peninsula were studied. The TEC from ROB-TEC within the seismogenic zone showed significant enhancement before the mainshock day. Ionospheric anomalies caused by solar activity and geomagnetic activity can be excluded, suggesting that earthquakes were the factor in ionospheric anomalies. Temporal TEC from ROB-TEC showed profound variations during the 3–10 days before the mainshock on relatively quiet storm days (Figure 4). Such evidence has already been provided by multiple analyses of temporal TEC over epicentral regions. They (Ulukavak (2020); Nico (2021); Nina (2020)) analyzed global earthquakes with Mw ≥ 6 and possible ionospheric TEC anomalies that occurred before earthquakes. The results of this study contribute to earthquake prediction in the short-term by providing empirical evidence [33–35]. Biagi revealed the INFREP Radio Network variations in connection with six earthquakes (Mw > 5.0) that occurred in the Balkan Peninsula and Adriatic Sea on 26 and 27 November 2019 [36]. All this information increases our knowledge about the origin of different pre-seismic signals associated with the above-mentioned catastrophic earthquakes.

In addition, with the increase in the number of seismogenic zones, the incidence and intensity of ionospheric disturbance has increased, further analyses are needed to quantify the variations with a more statistical and mathematical model and method to convincingly prove such phenomena. However, there are many theories regarding the existence of abnormal TEC from GPS measurements and some have provided enough evidence for the connection of the lithosphere–atmosphere–ionosphere coupling system. For example, Freund has proved the existence of stress-activated electric currents in rocks that sweep through the mineral grains under stress dislocations and cause the peroxy links to break, thus generating flow-down stress gradients, constituting an electric current with attendant magnetic field variations and EM emissions. Further experiments have shown

the generation of energetic charge barriers and their thermal conduction as fissure reactions in the seismogenic zone, which could be an intrinsic source of earthquake–ionospheric anomalies [37]. The main findings show that continuous earthquakes have repeatedly triggered the generation of high-energy charge barriers in the seismogenic belt and their heat conduction as crack reactions, which frequently affect the ionosphere.

On the other hand, CID observed by GPS stations (Figure 7) shows that with the increase in the number of earthquakes, the ionosphere over the seismogenic area is more obviously disturbed, which is consistent with the previous findings. For example, the farfield CID following the 2015 Mw 7.8 Nepal earthquake was detected by GPS-TEC [38] as the acoustic gravity wave induced by Rayleigh surface wave propagated into the ionosphere and caused the electronic and plasma density oscillation.

It is worth noting that after three strong earthquakes, the TEC suddenly decreased over the seismogenic area, forming a phenomenon similar to an ionospheric hole (Figure 6b). In further studies, more statistical and mathematical models and methods are required to prove such phenomena. There is evidence of a quantitative relationship between the TEC depression rate of the ionospheric cavity caused by the tsunami earthquake and the intensity of the tsunami earthquake [32,39]. However, more relevant observation data are needed to analyze the ionosphere at different heights in the region and the relevant meteorological data in the region are necessary for the generation of the ionospheric hole. In conclusion, the abovementioned results provide an opportunity to study the correlation between the ionosphere and earthquakes from multiple directions and angles for the earlier detection of specific pre-seismic anomalies related to earthquakes.

## **6. Conclusions**

In this study, we presented research on pre-seismic ionospheric anomalies and CID for a seismic swarm in the Balkan Peninsula in November 2019. The results show that a larger number of positive anomalies appeared in the seismogenic area on 23 November 2019 (3 days before the earthquake), and continuous positive anomaly disturbances appeared in the grid points in the spatially overlapping seismogenic zones. Moreover, there was a small range of negative anomalies in the seismogenic area 6–12 days before the earthquake, although the interval was relatively long. We found that CID could be clearly observed at the stations in the seismogenic area during each earthquake, and the degree of CID and disturbance increased as the IPPs approached the center of the earthquake area. After the three earthquakes, we also found a similar ionospheric hole phenomenon over the seismogenic area, and a large area where the VTEC value plummeted over the seismogenic area. The results showed that with the increase in the number of seismogenic zones, the incidence and intensity of the ionospheric disturbance increased.

These findings suggest that there is synergy between different surface, atmospheric, and ionospheric processes. However, further analysis is required to support the hypothesis regarding lithosphere and ionosphere coupling during the earthquake preparation period as mainshock responses, which will be left for a future work.

**Author Contributions:** Conceptualization, J.L. and L.L.; methodology, L.W. and J.L.; validation, C.R., L.H. (Liangke Huang) and X.T.; formal analysis, L.W.; writing—original draft preparation, L.W., L.Z. and D.Z.; writing—review and editing, D.Z., J.L., L.H. (Ling Huang) and H.H.; funding acquisition, J.L., D.Z., L.H. (Ling Huang) and L.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the Guangxi Natural Science Foundation of China (2020GXNS-FBA297145, 2020GXNSFBA159033), the Foundation of Guilin University of Technology (GUTQDJJ661- 6032), and the National Natural Science Foundation of China (42004025, 42064002).

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.
