Spatial analyses on pre-earthquake ionospheric anomalies and magnetic storms observed by China seismo-electromagnetic satellite in August 2018

The China Seismo‑Electromagnetic Satellite (CSES), with a sun‑synchronous orbit at 507 km altitude, was launched on 2 February 2018 to investigate pre‑earthquake ionospheric anomalies (PEIAs) and ionospheric space weather. The CSES probes manifest longitudinal features of four‑peak plasma density and three plasma depletions in the equato‑ rial/low‑latitudes as well as mid‑latitude troughs. CSES plasma and the total electron content (TEC) of the global ionosphere map (GIM) are used to study PEIAs associated with a destructive M7.0 earthquake and its followed M6.5 and M6.3/M6.9 earthquakes in Lombok, Indonesia, on 5, 17, and 19 August 2018, respectively, as well as to exam‑ ine ionospheric disturbances induced by an intense storm with the Dst index of − 175 nT on 26 August 2018. Anomalous increases (decreases) in the GIM TEC and CSES plasma density (temperature) frequently appear specifi‑ cally over the epicenter days 1–5 before the M7.0 earthquake and followed earthquakes, when the geomagnetic conditions of these PEIA periods are relatively quiet, Dst: − 37 to 19 nT. In contrast, TEC and CSES plasma parameter anomalies occur globally in the southern hemisphere during the storm days of 26–28 August 2018. The CSES ion velocity shows that the electric fields of PEIAs associated with the M7.0 earthquake are 0.21/0.06 mV/m eastward and 0.11/0.10 mV/m downward at post‑midnight/post‑noon on 1–3 August 2018, while the penetration electric fields during the storm periods of 26–28 August 2018 are 0.17/0.45 mV/m westward/downward at post‑midnight of 02:00 LT and 0.26/0.26 mV/m eastward/upward at post‑noon of 14:00 LT. Spatial analyses on CSES plasma discriminate PEIAs from global effects and locate the epicenter of possible forthcoming large earthquakes. CSES ion velocities are useful to derive PEIA‑ and storm‑related electric fields in the ionosphere.


Introduction
A sun-synchronous satellite, the China seismo-electromagnetic satellite (CSES, also named ZhangHeng-1), orbiting at 507 km altitude with ascending/descending node time of 02:00/14:00 LT (local time), was launched on 2 February 2018 to detect pre-earthquake ionospheric anomalies (PEIAs) and to study space weather effects in the ionosphere.Since satellite observations provide global and uniform coverage, they are ideally used to observe ionospheric weather and PEIAs (Pulinets and Boyarchuk 2004;Ouzounov et al. 2018).
For satellite remote sensing, Liu et al. (2001) pioneered applying the total electron content (TEC) derived from measurements of a local network of 13 ground-based global positioning system (GPS) receivers in Taiwan and found the GPS TEC anomalously decreases around the epicenter of days 1, 3, and 4 before the 21 September 1999 Mw 7.6 Chi-Chi earthquake.Following Liu et al. (2001), to detect PEIAs in a certain region, we can statistically examine and find the PEIA characteristics of increases or decreases (i.e., positive or negative polarity), appearance local time, duration, lead days, etc. of TEC anomalies appearing before large earthquakes over the region during a long-term period.After the characteristics are found, when TEC anomalies meet the characteristics, we might consider that temporal PEIAs have been detected in the region.Based on the found characteristics, scientists reported the temporal PEIA associated with large earthquakes over Taiwan (Liu et al. 2001(Liu et al. , 2004a(Liu et al. , b, 2010a;;Nishihashi et al. 2009), Indonesia (Liu et al. 2010b), China (Liu et al. 2009(Liu et al. , 2018a;;Chen et al. 2015), Japan (Kon et al. 2011;Liu et al. 2013Liu et al. , 2018b) ) and the USA (Su et al. 2013).However, these temporal PEIAs often suffer from global effects, such as solar disturbances, solar flares, magnetic storms, etc.
To discriminate local effects, such as PEIAs, from global ones, a spatial analysis of the TEC in the global ionosphere map (GIM) derived from the global navigation satellite system (GNSS: GPS, GLONASS, Galileo, or Beidou systems) is ideal to be utilized to find the latitude-longitude locations and/or distribution of TEC anomalies (cf.Sun et al. 2017).Since each map consists of 5183 (= 71 × 73) grid points (lattices), instead of examining the TEC over a single lattice of the monitoring region (i.e., the temporal analysis), the spatial analysis, in fact, allows us simultaneously examine the latitudelongitude distribution of TEC anomalies over the global 5183 lattices.When anomalies frequently appear over a small specific area, it shows that a local effect, such as an earthquake-related anomaly, has been observed.Moreover, if the anomalies detected by the spatial analysis also meet the PEIA characteristics at the region, we might then consider the PEIA being observed.Several articles report spatial PEIAs of GIM TEC associated with the 16 October 1999 Mw7.1 Hector Mine earthquake (Su et al. 2013), the 2004 M9.3 Sumatra-Andaman earthquake (Liu et al. 2010b), the 2008 M8.0 Wenchuan earthquake (Liu et al. 2009), the 2010 M7.0 Haiti earthquake (Liu et al. 2011), and the 11 March 2011 M9.0 Tohoku earthquake (Liu et al. 2018b).Since the spatial analysis conducts a global search on the location and frequency of anomalies, the spirit is similar to simultaneously examining and/or testing anomalies in 5183 possible earthquake locations.By contrast, when anomalies worldwide occur, it would result from global effects, for instance, magnetic storms.Therefore, the spatial analysis can be used to separate local effects from global ones, verify the temporal PEIA being detected, and locate possible forthcoming large earthquakes.
For a satellite in situ measurements, DEMETER (Detection of Electro-Magnetic Emissions Transmitted from Earthquake Regions) might be the first satellite designed specifically to in situ detect PEIAs, which is a micro-satellite that investigates the region within 65° geomagnetic latitude at a 680 km altitude, 98.3° inclination orbit during June 2004-December 2010.Many research articles report PEIAs in the electron density, electron temperature, ion density, and ion temperature (Ho et al. 2013a(Ho et al. , b, 2018;;Shen et al. 2015Shen et al. , 2017;;Liu et al. 2015Liu et al. , 2016;;Tao et al. 2017;Yan et al. 2017), as well as amplitude changes of electro-magnetic emissions (http:// smsc.cnes.fr/ DEMET ER/A_ publi catio ns.htm, cf.Nemec et al. 2008).In addition to the scalar quantities, such as electron density, electron temperature, ion density, ion temperature, etc., Liu and Chao (2017), for the first time, use vector quantities of the ion velocity probed by ROC-SAT to estimate seismo-generated electric field appearing 2 days before the 31 March 2002 M6.8 Earthquake in Taiwan.Based on Liu and Chao (2017), Ho et al. (2018) analyzed DEMETER observations before 49 M ≥ 6.5 earthquakes in 2010 and showed that the perpendicular ion velocity and associated electric field could be anomalously changed before the earthquakes.

Observations and data analyses
The CSES satellite carries eight instruments, including a search-coil magnetometer, electric field detector, high precision magnetometer, GNSS occultation receiver, plasma analyzer package (PAP), Langmuir probe (LAP), high energetic particle package and detector, and tri-band beacon (Lin et al. 2018;Shen et al. 2018;Yan et al. 2018;Yang et al. 2020).Figure 1 illustrates global smoothed median maps of the CSES electron/ion density (Ne/Ni) and temperature (Te/Ti) observations at 14:00 and 02:00 LT in August 2018.Regarding ionospheric weather, Fig. 1a-d displays along the magnetic equator that manifest longitudinal features of four Ne/Ni peaks (wavenumber four, WN4) (Sagawa et al. 2005;Immel et al. 2006;Lin et al. 2007;Kil et al. 2015) at post-noon of 14:00 LT and three plasma depletions (plasma depletion bay, PDB) (Chang et al. 2020(Chang et al. , 2024) ) at post-midnight of 02:00 LT.The WN4 occurs centering at near − 160E, − 90E, 20E, and 110E longitude, while the PDB appears at around -30°E, 60°E, and 150°E longitudes.Owing to the cooling through Coulomb collisions (Kakinami et al. 2011), Fig. 1e-h further depicts along the magnetic equator that longitudinal variations of Te/Ti are generally in opposite polarities with those of Ne/Ni.For WN4 (PDB), the four maxima (three minima) of Ne/Ni are in coincidence with the four minima (three maxima) of Te/Ti.Moreover, Fig. 1c, d, g, and h reveals that post-midnight Ne/Ni (Te/Ti) significantly reduce (suddenly enhance) within 45°-65° geographic latitude, which shows features of the mid-latitude trough (cf., Lee et al. 2011;Chen et al. 2018).
Table 1 displays the occurrence time and location of the M7.0 earthquake on 5 August and its three followed M6.s earthquakes on 17 and 19 August 2018.Figure 2a from top to bottom displays the solar radio flux at 10.7 cm of F10.7, and the Kp and Dst magnetic indices, which show the magnetic condition being generally quiet during 1-25 August 2018.The Dst index displays an intense storm that reached the maximum depression of -175 nT on 26 August 2018.The CODE (Center for Orbit Determination in Europe) GIM of the TEC maps, with 1-h time resolution (24-timepoints, at the top of the hour daily) covering ± 87.5°N latitude and ± 180°E longitude ranges with spatial resolutions of 2.5° and 5°, respectively, is ideal to be utilized for detection of TEC anomalies.Each map consists of 5183 (= 71 × 73) grid points (lattices).TECs along the epicenter longitude and over the epicenter of the M7.0 earthquake were extracted from a sequence of the GIM images in August of 2018 to study anomalies.On 26 August 2018, Fig. 2b shows the latitude-time-TEC (LTT) plots along the epicenter latitude and displays that the EIA (equatorial ionization anomaly) crests move poleward and yield the greatest TEC values, which suggests that the daily dynamo electric field being superimposed with the prompt penetration electric field (Kelley 2009;Liu et al. 2013) significantly increases.
We further detected abnormal signals of TEC variations over the epicenter by applying a quartile-based process in August 2018.For each timepoint of the day, we compute the moving median M 30 days before the observation day as well as find the deviation between the observed values and the computed median.To provide information about the deviation, we also calculate the first (or lower) and the third (or upper) quartiles, denoted by LQ and UQ, respectively.To have a stringent criterion, we set the threshold k = 1.5 to construct the lower bound, LB = M −k ( M−LQ), and the upper bound, UB = M + k (UQ− M ).Therefore, the probabil- ity of a new TEC in the interval (LB, UB) is approximately 65% (Klotz & Johnson 1983).The median together with the associated LB and UB, then provides references for detecting abnormal signals in TEC variations on the observation day.When an observed TEC is not in the associated (LB, UB), we declare a positive or negative abnormal TEC signal.Figure 2c depicts that pronounced positive TEC anomalies (red shaded areas) appeared day 2, 4 prior to the 5 August M7.0 earthquake on 1, 3 August and days 1-5 before the 17 August M6.5 and the 19 August M6.3/M6.9earthquakes on 12, 14-15, 17-18 August (also see Table 1).
The positive TEC anomalies appearing on 12, 14, 15, 17, and 18 August 2018 seems to be post-seismic signatures on days 7-13 after the M7.0 earthquake, but possibly to be PEIAs occurring days 1-5 before the M6.5/ M6.3/M6.9earthquakes.In general, the sequency of possible PEIAs tends to lead that of the earthquakes by about 1-5 days (Table 1).The most prominent positive TEC anomalies occurred during the storm period of 26-28 August.Note that almost no prominent negative anomalies (black shaded areas) have been detected in August 2018.If the above-observed positive TEC anomalies several days before the four large (M7.0/M6.5/M6.3/M6.9)earthquakes in August 2018 were indeed PEIAs, we could expect those past large earthquakes in the same area were often preceded by similar anomalies with similar lead times.We first find the characteristics of PEIAs of GIM TEC associated with 123 M ≥ 6.0 earthquakes occurring over the Lombok area (6.7°N-23°S by 101.4-131.4°E)during the 11-year period from 2007 to 2017 (Fig. 3).We then apply the z test (Chen et al. 2015;Liu et al. 2018bLiu et al. , 2023) ) to examine the GIM TEC associated with the 123 earthquakes to find the characteristics of M ≥ 6.0 earthquakes over the Lombok area.The z value is then used to find if at certain timeperiods, the observed proportion of ) and peaks 74 fsu on 24 August.The Dst index displays that it is relative quiet, especially 1-10 August, and there is an intense magnetic storm with a maximum depression of − 175 nT (Kp 7 + ) on 26 August.The LTT plots reveal that the EIA (equatorial ionization crest) becomes most intense and moves very poleward on the storm day of 26 August.The most significant positive TEC are observed before those earthquakes on 1, 3, 12, 14, 15, 17, and 18 August 2018 earthquake-related anomalies differs significantly from the background (i.e., normal time) proportion of anomalies.For each timepoint of the day, let π be the observed proportion of earthquake-related anomalies, which is the number of positive or negative anomalies derived by the number of the M ≥ 6.0 earthquakes (123 earthquakes), and π 0 be the background proportion of anomalies in the 11 years of 2007-2017 period (4018 days), which is the number of positive or negative anomalies derived by the number of days in the 11-year period (4018 days).The z value is then given by where n = 123 is the number of earthquakes.If z > 1.96 , we claim, at a significant level of 0.025, that π > π 0 .The z test is conducted for positive and negative anomalies separately.Figure 4a displays that the positive anomalies (red contours) appear frequently between day 8 before and day 4 after, especially on day 4 before and after, the 123 M ≥ 6.0 earthquakes, which generally agrees with the positive anomalies (red shaded areas) appear day 2, 4 (1-5) before the M7.0 (M6.5/M6.3/M6.9)earthquake on 1, 3 (12, 14-15, 17-18) August 2018 shown in Fig. 2c.On the other hand, Fig. 4b reveals negative anomalies (blue contours) rarely occur between day 13 before and day 14 after the 123 M ≥ 6.0 earthquakes, which also agrees with that almost no obvious negative anomalies appear several days before and after the four earthquakes shown in Fig. 2c.In Fig. 4a, c, a positive TEC anomaly zone with statistical significance at 05:00-09:00 UT on day 4 before the 123 M ≥ 6.0 earthquakes indicates that the PEIA characteristics of M ≥ 6.0 earthquakes over the Lombok area are the positive polarity in GIM TEC, the local time period of 05:00-09:00 UT, and the lead time of about 4 days.
To verify if the positive TEC anomalies (or PEIAs) are a candidate for a reliable earthquake precursor, we treat the positive TEC anomalies as earthquake alarms at 05:00-09:00 UT day 4 before the M ≥ 6.0 earthquakes (Zone A in Fig. 4c), and construct the receiver operating characteristic (ROC) to evaluate the reliability of the earthquake alarming.Based on the PEIA characteristic for k = 1.5, when positive TEC anomalies appear more than onethird of the period of 05:00-0900 UT, we can issue an alarm for an earthquake with the magnitude of M ≥ 6.0 occurring in the following 4 days.Note that k = 1.5 is simply an imperial threshold in the previous studies (Chen et al. 2015;Liu et al. 2018aLiu et al. , b, 2023)).In fact, for small k values, many alarming days are issued, and hence more false alarms are obtained.However, for large k values, limited alarming days are issued and both the false alarm and successful rates are drastically reduced.Hence, to test the preference of the positive anomalies, we consider constructing the ROC curve for various k values.For each k value, we examine four different conditions, an alarm day being followed by earthquakes or no earthquake, and a non-alarm day being followed by earthquakes or no earthquake within a certain lead day period.Let TP(k) and FN(k) stand for the numbers of earthquake days with and without being led by alarm days, respectively.Moreover, denote FP(k) and TN(k) to be numbers of non-earthquake days with and without being alarmed, respectively.Then, the true positive rate TPR(k) and false positive rate FPR(k) as given by and where TPR(k) is the probability that an earthquake is successfully alarmed, and FPR(k) is the probability to make a false alarm.Hence, the ROC curve with FPR(k) as the x-axis and TPR(k) as the y-axis can be constructed.The area under the ROC curve (AUC) is further used for assessing the effectiveness of the PEIAs (Hanley and McNeil 1982).Note that when TPR(k) = FPR(k) for all k, an equal chance to alarm earthquake day and nonearthquake day, we have AUC = 0.5.Therefore, a reliable precursor should have AUC > 0.5.Note that the value of k varies from 0 to 10 by increasing 0.1, and therefore, there are 101 k values under investigation.Figure 5 shows that the ROC curves (color curves) of the positive TEC anomalies, a TPR is generally greater than the associated FPR, and the AUC is greater than 0.5 in the zone.The curve is beyond the associated 95% confidence; the AUC value of 0.54 is greater than 0.5; and the p value is zero, which confirms that the GIM TEC significant increases (positive anomalies) at 05:00-09:00 UT day 4 before meets the characteristic in detecting the 123 M ≥ 6.0 earthquakes in the Lombok area (Fig. 4).Therefore, the positive TEC anomalies appearing on 1-3 August 2018 (4-2 days before) can be considered the temporal PEIAs of the 5 August 2018 M7.0 Lombok earthquake.For the general/imperial index: k = 1.5, TP = 51, FP = 1133, FN = 72, TN = 2743, alarm days: 1185.Based on the above, the success rate, the number of earthquakes preceded by PEIAs divided by that of the total earthquakes, is 0.41 (= 51/123); the alarm rate, the number of alarm days divided by that of the total study days, is 0.29 (= 1185/4018); and the probability gain is of about 1.41 (= 0.41/0.29).For Youden index: K = 1.8, TP 49, FP = 944, FN = 74, TN = 2932, alarm days: 993; the probability gain is of about 1.60 (= (49/123)/(993/4018)).These show that PEIAs of positive (increase) TEC anomalies is a reliable earthquake precursor, and the empirical k = 1.5 is a good choice to have the maximum R score.
Figures 4 and 5 show that PEIAs of positive TEC anomalies frequently appear days 1-4, especially day 4, before the M ≥ 6.0 earthquakes, while the pre-earthquake TEC anomalies for the M7.0 + M6.s earthquakes seem to have p-value is "zero" occurred 1-5 days before the four earthquakes (Fig. 2c and Table 1).Hence, in the subsequent analyses of possible PEIAs of these four earthquakes, we assume a lead time range of 1-5 days.This accords with earlier studies of PEIAs (Song et al. 2020).Note this is a conservative choice considering that PEIAs tend to be more conspicuous as the time of earthquake approaches (e.g., Fig. 4 of Li et al. 2020).
To further augment the likelihood that the anomalies which preceded the M7.0 earthquake and its three M6.s earthquakes were actually PEIAs, as well as signatures of the magnetic storm, spatial analysis on the global distribution of positive TEC anomalies is conducted.Similar to the procedure in Fig. 2c, but for each lattice of GIM, we identify positive TEC anomalies and compute the proportion of the positive anomaly during the observation period, which is the number of identified positive anomalies being divided by the total spatial observation points (or lattices) of 5183.Thus, we apply the spatial analysis to globally find the frequency of positive anomalies associated with 5183 possible earthquake locations and compute the global distribution of the anomalies during PEIA periods of the M7.0 earthquake and three M6.s earthquakes.Figure 6a displays the proportions of the positive anomalies of the globe during the temporal PEIAs of the M7.0 earthquake on 1-3 August 2018 (72 timepoints in total), while Fig. 6b highlights that the large proportions specifically appear around the epicenter.Figure 6b depicts that the positive anomalies specifically appear most frequently over the epicenter (a rectangular of 25°S-5°N by 110°E-150°E) area, 0.1% (= 5/5183) of GIM, for more than 65% (= 47/72) of the PEIA time period, which shows that PEIAs of the M7.0 earthquake have been observed on 1 and 3 August 2018.Similar to Fig. 6a, b, for comparisons with the CSES observations, we compute proportions of the positive TEC anomalies of the globe associated with the three M6.s earthquakes of M6.5, M6.3, and M6.9 at 02:00 and 14:00LT on 12, 14, 15, 17, and 18 August 2018 (Fig. 6c), as well as associated with the M7.0 earthquake and its M6.s earthquakes (Fig. 6d). Figure 6c, d illustrates that the highest frequencies (60% = 6/10 and 64% = 9/14) of positive TEC anomalies are limited to the narrow (0.02% = 1/5183 and 0.06% = 3/5183) regions near the epicenters of the M6.s earthquakes and the M7.0 + M6.s earthquakes.Figure 6b-d illustrates that the top proportion of the anomaly duration associated with the M6.s earthquakes of 60% is smaller than that associated with the M7.0 earthquake of 65% and M7.0 + M6.s earthquakes of 64%, which shows that the greater earthquake appears more frequently in the PEIA period.
Owing to the quasi neutrality of plasma, both Ne and Ni are approximately proportional to TEC.Therefore, positive CSES Ne and Ni anomalies most likely significantly increase on the seven PEIA days of 1, 3, 12, 14, 15, 17, and 18 August 2018.On the other hand, due to the cooling through Coulomb collisions, CSES Te and Ti very likely decrease anomalously on the 7 days.Similar to Fig. 6e, we examine the spatial distribution of anomalous increases in Ne and Ni as well as abnormal decreases in Te and Ti at 14:00 and 02:00 LT on the seven PEIA days in taking the rest of 21 (= 31-7-3) days in August as a reference, except the storm days of 26-28 August 2018.Figure 7a, b shows during the PEIA days that the post-noon Ne and Ni tend to anomalously increase and specifically appear over the epicenters.Figure 7c, d depicts that the post-midnight Ne and Ni incline to abnormally increase mainly around the equatorial PDB longitudes and the mid-latitude trough zones.Note that the post-midnight Ni also increases abnormally over the epicenters (Fig. 7d).By contrast, Fig. 7e, f shows that the post-noon Te and Ti anomalously decrease specifically over the epicenters.Figure 7g, h illustrates that the post-midnight Te and Ti abnormally decrease around the mid-latitude trough zones.In fact, the post-midnight Ti also abnormally decreases around the equatorial PDB longitudes (Fig. 7h).In general, the CSES Ne, Ni, Te and Ti at the post-noon seem to have a better chance than those at the post-noon detecting PEIAs.Figure 7a-h illustrates that five anomaly types of the post-noon Ne increases, the post-noon Ni increases, the post-midnight Ni increases, the post-noon Te decreases, and the post-noon Ti decreases occur over the epicenters during the seven PEIA days.We father find the global distributions of anomaly types of CSES plasma parameters associated with the M7.0 earthquake, the M6.s earthquakes, and the M7.0 + M6.s earthquakes.It can be seen that the top count of types associated with the M7.0 earthquake (Fig. 7i) and the M6.s earthquakes (Fig. 7j) scatter around the epicenter area, while those associated with the M7.0 + M6.s earthquakes (Fig. 7k) appear in the narrow region near the epicenters.Agreements of these anomalies frequently appearing over the epicenters in Figs.6b-d and 7i-k show that PEIAs can be observed by GIM TEC and CSES plasma parameters, while the PEIA is associated with the M7.0 earthquake on 1-3 August is the most intense one.The spatial analyses in Figs.6d and 7k show that the GIM TEC and CSES plasma parameter anomalies frequently occur in the narrow region near the epicenters day 4 before the M7.0 + M6.s earthquakes.
Figure 6e illustrates the proportions of the positive anomalies of the globe during the storm period on 26-28 August 2018, while Fig. 6f reveals that the large proportions (more than 75%) occur in the southern hemisphere, especially along the magnetic equator ± 15°, 20-40°S magnetic latitude, and 45-75°S magnetic latitude.Positive anomalies frequently appear in most of the southern hemisphere, which confirms the global effect that a positive ionospheric storm occurs during 26-28 August 2018.Similar to Fig. 7, we examine the spatial distribution of anomalous increases in Ne and Ni as well as abnormal decreases in Te and Ti during the storm period of 26-28 August 2018.Figure 8 displays that the most frequent location of storm-related anomalies, which can be obtained by anomalies in the post-noon/post-midnight Ne, Ni, Te, and Ti (Fig. 8a-h), is mainly in the southern hemisphere.It is worth noting that Fig. 8i which is the distribution of three or more anomaly type occurrences and Fig. 6f are in excellent agreement.
Based on the above good agreement between Figs. 6i and 8f, we then focus on CSES data within the rectangular area (25°S-5°N by 110°E-150°E) in Figs. 6, 7, and 8 to have a better understanding of the causal mechanisms Since ionospheric data inhabit a right-skewed and/or heavy-tailed distribution, it is suitable to apply a median base analysis.The box-and-whisker procedure (Wilcox 2010), as a median base analysis, has the advantage of visually observing the significant difference among multi-data sets simultaneously.Therefore, we apply boxand-whisker (box) plots to statistically investigate the CSES V D , V E , and V N data within the rectangular area at the post-midnight and the post-noon on the observation and the reference days (i.e., all of the days in August other than 1-3, 12, 14, 15, 17, 18, and 26-28) for the PEIA period of the M7.0 earthquake on 1-3 August and the M7.0 earthquake plus the three followed earthquakes on 1-3, 12, 14, 15, 17, and 18 August, as well as the storm period of 26-28 August.The ends of a box denote LQ and UQ.The horizontal line within the box denotes the median.If two boxes do not overlap with each other, we claim that there is a dramatic difference between the two boxes.However, when one shorter box with a median larger than the upper quartile or smaller than the lower quartiles of the other longer box, the two boxes are still declared to be different.Figure 10 displays the boxes of the associated references (gray boxes) and the observations (black boxes) of the M7.0 earthquake (mainshock), M6.s earthquakes, and M7.0 + M6.s  (a-h) for the lattice being equal to and/or greater than 3 counts.The criteria for anomalies are the same as Fig. 7 (Te/Ti < 25% and Ne.Ni > 50%).It can be seen that the most prominent counts mainly occur in the southern hemisphere earthquakes in CSES plasma drifts (ion velocities).Based on the dynamo theory (cf.Kelley 2009), the ionospheric electric field E = − V × B, can be estimated, where V is the ion velocity and the B is the Earth's magnetic field.From the IGRF (international geomagnetic reference field) model, we find the B field at the satellite orbit height of 507 km altitude over the epicenter is 4.7 × 10 -5 T with a magnetic dip of -39.44 degrees and the declination of 1.64 degrees over the M7.0 earthquake area.Owing to the relatively low magnetic dip of the earthquake location and the high satellite speed in the meridional direction, we focus on the eastward and downward ion drifts and derive the vertical and zonal electric fields, respectively.By giving the median value of the velocities of the two boxes, we can calculate the electric fields on the observation and reference days, as well as their difference, subtracting the former from the latter.Figure 10 reveals that velocity differences between PEIA days and reference generally are insignificant, except that the post-midnight V D = − 11.2 m/s of the M7.0 earthquake PEIA days on 1-3 August 2018 significantly differs from the reference of V D = − 3.8 m/s.Then, the PEIA-related electric field Fig. 9 GIM TEC as well as CSES Ne, Ni, Te, Ti, V D , V E , and V N in the epicenter rectangular (25°S−5°N by 110°E−150°E) area.Data at 02:00 LT (left panels) and at 14:00 LT (right panels).The red and blue dashed lines denote the PEIA period of 1-3 August 2018 and the storm of 26-28 August 2018, respectively.The star symbols denote the M7.0, M6.5, M6.3, M6.9 earthquake, respectively (See figure on next page.)Fig. 10 Box plots of the associated reference (gray box in each panel) and the observation (black box in each panel) of the downward velocity (top row panels) and eastward velocity (second row panels) at 0200 LT (left column panels) and 14:00 LT (right column panels) during the M7.0 earthquake (first right-side black boxes), M6.s earthquakes (second right-side black boxes), and M7.0 + M6.s earthquakes (very right-side black boxes) PEIA periods.Box plots of the storm period are given in third row and bottom panels.The horizontal line within the box denotes the median.The ends of the box are the first quartile (25% of the data set, Q1) and third quartile (75% of the data set, Q3), where the first (third) quartile is the middle value between the smallest (highest) and the median of the data set.The difference between the Q1 and Q3 is called the inter-quartile range (IQR).If a value lower than Q1−1.5IQR(lower red dot) and/or greater than Q3 + 1.5IQR (upper red dot), it is declared the outlier (black dot).The horizontal lines extending out from the box are the minimum and maximum value, which the minimum (maximum) is the smallest (largest) value within the range of outlier.The vertical lines out from the box to the minimum (maximum) are called the lower (upper) whiskers.CSES Ne (Ni), Te (Ti), and ion velocities (V D , V E , and V N ) have been underestimated by a factor of 10, underestimated by about 1300 °K, and overestimated by a factor of 20, respectively.Note that the calibrated and raw data are, respectively, giving in left-side and right-side axes ΔE = 0.21 mV/m eastward can be obtained by inserting the velocity difference of − 7.4 (= − 11.2 + 3.8) m/s into E = − V × B. Although differences in the rest velocities are insignificant, we find the direction of the PEIA-related electric fields associated with the M7.0 earthquake PEIAs being 0.11 mV/m downward at post-midnight and 0.06 mV/m eastward and 0.10 mV/m downward at postnoon.The PEIA-related electric fields at post-midnight/ post-noon tend to be eastward and downward.In contrast, the velocity differences are generally significant for the storm days, except for the post-noon V E .The storm electric fields are 0.17 mV/m westward and 0.45 mV/m downward at post-midnight, and 0.26 mV/m eastward and 0.26 mV/m upward at post-noon.

Conclusion and discussion
Figures 1 and 7 reveal the ionospheric weather features of the CSES Ne, Ni, Te, and Ti in WN4s, PDBs, and midlatitude troughs.Figures 2 and 6, 7, 8 show that PEIA and storm signatures can be detected by both GIM TEC and CSES plasma quantities.These indicate that CSES well meets its mission goals of detecting seismo-electromagnetic signals and observing space weather.The spatial analyses of GIM TEC and CSES plasma quantities allow us to discriminate the local anomalies of PEIAs (Figs. 6 and 7) from the global anomalies of WN4s, PDBs, mid-latitude troughs, magnetic storms (Figs.6f and 8i), etc., which can be used to not only verify the temporal PEIAs being detected but also locate possible forthcoming earthquakes.In addition to the statistical evidence of PEIA existence (Figs. 2c,4,5,6 and 7), the CSES velocities allow us to have a better understanding of the causal/ physical mechanisms of the detected PEIAs.Velocity differences between the M6.s or M7.0 + M6.s PEIA days and the reference is insignificant, and however, V D between the M7.0 PEIA days and the reference is significantly different (Fig. 10).Meanwhile, Ho et al. (2018) analyze both ion density and ion velocity observed by DEMETER before 49 earthquakes with M ≥ 6.5 in 2010, and find that both ion density and ion velocity increase anomalously before the earthquakes.They further report the differences before and after earthquakes are more significant with larger magnitude, and the density variations have a positive correlation with the perpendicular velocity, which is related to the perpendicular dynamo electric field.Results of Fig. 10 and Ho et al. (2018) show that the velocity difference is proportional to the earthquake magnitude, and the greater earthquake is preceded by the larger (smaller) PEIAs of GPS TEC, Ni, and Ne (Ti and Te).Based on the dynamo theory, we find that PEIArelated electric fields associated with the M7.0 earthquake PEIAs are 0.21 mV/m eastward and 0.11 mV/m downward at post-midnight and 0.06 mV/m eastward and 0.10 mV/m downward at post-noon.Note that the PEIA-related electric fields remain in the same eastward and downward directions at both post-midnight and post-noon.It is most likely that due to the eastward PEIA-related electric fields, E × B upward plasma drifts uplift the ionosphere, which in turn results in the Ne and Ne increases as well as the Te and Ti decreases during the M7.0 earthquake PEIA days of 1-3 August 2018.It might be that the downward PEIA-related electric fields cause the eastward plasma drifts, which further result in the positive anomalies of GIM TEC occurring in the eastward side of the epicenter shown in Fig. 6a-d.
To explain the penetration electric field from high-to low-latitudes during storm periods, Kelley (2009) proposed that the increasing activity currents (Region 1 > Region 2) yield an eastward perturbation electric field on the dayside and a westward perturbation one at night.The storm electric fields of 0.17 mV/m westward and 0.45 mV/m downward at post-midnight, as well as those of 0.26 mV/m eastward and 0.26 mV/m upward at post-noon (Fig. 10) indicate that the increasing activity currents play an important role during the storm period of 26-28 August 2018.Notice that the storm perturbation electric fields yield opposite directions in westward/ downward at post-midnight and eastward/upward at post-noon.The same and opposite directions at postmidnight and at post-noon indicate that the generations of the PEIA-and storm-related electric fields are very different.The former will be generated around the epicenter area (local scale) during the PEIA period, while the latter is caused by the magnetosphere-ionosphere currents (global scale) during the storm period (for a complete overview see, Piersanti et al. (2020).Yang et al. (2020) examined time series data recorded by the magnetometer, Langmuir probe, plasma analyzer, electric field detector, and GNSS occultation receiver onboard the CSES satellite during the 25 August 2018 magnetic storm.They found that the storm possibly induced electric field penetration, causing the positive ionospheric storm, excited significant ELF/VLF waves, and enhanced the energetic electron flux.In contrast, we construct global maps of CSES Ni/Ne and Ti/Te together with GIM observations, which observe the global distribution of WN4, PDB, and mid-latitude trough features and storm signatures (Figs. 1,6,and 8).Figures 6f and 8 show that positive anomalies pronouncedly frequently appear in the southern hemisphere, which confirms the global effect that the positive ionospheric storm occurs during 26-28 August 2018.This positive ionospheric storm agrees well with Lissa et al. (2020) and Moro et al. (2021) that the neutral composition in [O]/[N2] observed by TIMED/GUVI remarkably enhanced in the southern hemisphere during the storm periods, especially on 26 August.The positive storm effect during the main and recovery phases is mainly attributed to the eastward prompt penetration electric fields in addition to a strongly enhanced ratio of thermosphere neutral composition.Here, we further apply the dynamo theory on the CSES V D and V E (Fig. 9) evaluating the strength of prompt penetration electric fields (Fig. 10).Song et al. (2020) applied the superposed epoch analyses method on GIM TEC associated with 35 M ≥ 5.8, depth ≤ 50 km earthquakes in Indonesia during 2007-2017 and find that the PEIA characteristics are significant positive anomalies appearing 3-5 and 7 days before the earthquakes.This generally agrees with the characteristics of PEIAs that positive anomalies in GIM TEC appear on day 4 before the 123 M ≥ 6.0 earthquakes in the Lombok area from 2007 to 2017 (Fig. 4).For event studies, Song et al. (2020) investigated PEIAs of four successive shallow-focus earthquakes (M6.4 on 28 July 2018, M6.8, 5.9, 6.9 on 05, 09, and 19 August 2018) in Indonesia by CSES Ne and Te as well as GIM TEC.They found remarkable CSES Ne enhancements 1, 5, 2 and 5 days before the four earthquakes.Liu et al. (2006) investigates the relationship between 184 M ≥ 5.0 earthquakes and PEIAs of the ionospheric F2-peak plasma frequency in the Taiwan area during 1994-1999.They find that the earthquakes and PEIAs appear sequency-by-sequency in tandem; in certain time windows with fewer earthquakes, the anomalies seem to appear less frequently; and the PEIAs lead the earthquakes by 1-5 days.These show that a PEIA sequency is tightly followed by an earthquake sequency, which could often result in one earthquake's PEIAs being the other earthquake's post-seismic signatures.Figure 2c depicts the sequency of positive TEC anomalies (or PEIAs) on 1, 3,12,14,15,17,18 August 2018 and the sequency of the 5 August M7.0,18 August M6.5, 19 August M6.3/M6.9earthquakes.Although they seem to be post-seismic signatures of the 5 August M7.0 earthquake, the positive TEC anomalies on 12, 14, 15, 17, and 18 August are in fact the PEIAs associated with the M6.5/M6.3/M6.9earthquakes (Table 1, Figs. 6c and 7j).The positive TEC anomalies (or PEIAs) sequency leads the M/7/M6.5/M6.3/M6.9earthquake sequency by about 1-5 days, which agrees with Song et al. (2020).
For spatial PEIA analyses, Fig. 6c, e depict that positive GIM TEC anomalies specifically appear around the epicenter 2-4 days before the 5 August 2018 M7.0 earthquake and 1-5 days before the three M6.s earthquakes, respectively, while Song et al. (2020) report that positive GIM TEC anomalies specifically appear around the epicenter 5 days before the 19 August 2018 M6.9 earthquake (one of the M6.s earthquakes in this study).Anomalies frequently appearing specifically around the epicenter a few days before large earthquakes in the two studies indicate that the electromagnetic environment around the epicenter of forthcoming large earthquakes might have been significantly modified.Song et al. (2020) examined the CSES Ne within 60°N-60°S by 0-180°E and found that the Ne abnormally increases 1-5 days before the four earthquakes, while Fig. 7k shows that the anomalies at post-noon Ne, post-noon Ni, post-midnight Ni, postnoon Te, and post-noon Ti appear around the epicenter 1-5 days before the M7.0 + M6.s earthquakes.Note that the location of the positive anomalies in the GIM TEC (Fig. 6a-d) and that of the anomalies in the CSES plasma parameters (Fig. 7k) are in good agreement.Li et al. (2020) examined the primary joint statistical seismic influence on ionospheric parameters recorded by the CSES and DEMETER satellites and showed that the CSES could effectively register ionospheric perturbations due to strong earthquakes as the DEMETER satellite does.They find that the detection rate of pre-earthquake Ni and Ne perturbations increases as the earthquake magnitude increase and as the focal depth decreases, which generally agrees with statistical results in studying pre-earthquake anomalies of the ionospheric F 2 -peak plasma frequency, foF 2 (Liu et al. 2006) and TEC (Liu et al. 2009(Liu et al. , 2018b;;Le et al. 2011).In fact, the maximum proportion of 65% of the M7.0 earthquake (Fig. 6b) being greater than that of 60% M6.0 s earthquakes (Fig. 6c) and 64% of the M7.0 + M6.s earthquakes (Fig. 6d) also show the larger earthquake, the better chance for the earthquake to be recognized by the anomalies.Liu et al. (2015) reported that the intersections of the global distribution of anomalies in the nighttime Ne, the nighttime Ni, and the daytime Ti of DEMETER frequently appear specifically over the epicenter 1-6 days before the 2008 M8.0 Wenchuan earthquake.Similarly, Fig. 7 shows that the global distribution of the top count of anomalies in CSES post-noon/post-midnight Ne, Ni, Te, and Ti also frequently appear in the narrow region near the epicenters 1-5 days before the M7.0 + M6.0 s earthquakes.Anomalies specifically appear over the epicenter a few days before large earthquakes in Liu et al. (2015) and Fig. 7 again shows that both CSES and DEM-ETER can be used to detect PEIAs and the spatial analysis locate the epicenter of possible forthcoming large earthquakes.
In summary, the GIM TEC and the CSES ion density are used to study PEIAs associated with the M7.0, M6.5, and M6.3/M6.9earthquakes in Lombok, Indonesia on 5, 17, and 19 August 2018, respectively, as well as to examine ionospheric disturbances induced by the intense storm on 26 August 2018.The significant positive TEC anomalies associated with the M7.0 and M6.s earthquakes appeared on 7 days of 1, 3, 12, 14, 15, 17, and 18 August 2018.The spatial analyses show that the positive TEC anomalies frequently occur specifically near the epicenters, which indicates that the PEIAs of the M7.0 and M6.s earthquakes have been observed.The CSES Ne, Ni, Te, and Ti are used to detect PEIAs and to observe ionospheric weather, such as WN4, PDB, and mid-latitude trough, magnetic storms, etc.The spatial analyses allow us to discriminate the PEIA from global signatures of ionospheric storms, WN4, PDB, mid-latitude trough, etc.Similar tendencies in concurrent and co-located measurements of the GIM TECs and the CSES Ne and Ni indicate that the two observations can be used to threedimensionally detect PEIAs and examine ionospheric storm signatures.The CSES ion velocity has a better understanding of the causal mechanisms of the PEIAand storm-related electric fields.The PEIA-related electric fields of 0.21 mV/m and 0.06 mV/m eastward appear, respectively, at post-midnight and post-noon on day 4 before the M7.0 earthquake, while the prompt penetration electric fields of 0.26 mV/m eastward and 0.26 mV/m downward at post-noon as well as 0.17 mV/m westward and 0.45 mV/m downward at post-midnight are estimated.In conclusion, results demonstrate that spatial analyses on CSES plasma parameters can be employed to detect PEIAs of large earthquakes, and observe ionospheric space weather of the globe.

Fig. 1
Fig. 1 Global observations of CSES in a-d Ne/Ni and e-h Te/Ti at 14:00 LT and 02:00 LT in August 2018.The red star denotes the epicenter of the 5 August 2018 M7.0 earthquake.CSES data were first grided into 5183 lattices with a lattice size of 2.5° by 5° in latitude by longitude.The median value of each lattice is computed, and the medians are smoothed for the global map by a window of 5-lattice by 5-lattice sliding by 1-lattice in the latitudinal and longitudinal directions

Fig. 2
Fig. 2 Solar radio flux, magnetic condition, and ionospheric TEC variations in August 2018.a From top to bottom are the solar radio flux at 10.7 cm of F10.7 and the magnetic indices of Kp and Dst.b Latitude-time-TEC (LTT) plots along the epicenter longitude.c Time series of the GPS over the epicenter in August 2018.The 5 August 2018 M7.0 earthquake and its followed M6.s earthquakes of 17 August M6.5 and 19 August M6.3/M6.9 are denoted by the vertical red solid and dash lines, respectively.The red curve is the observed GIM TEC.The gray and black curves demote the associated median ( M ), upper bound (UB), and lower bound (LB), respectively.The red and blue shaded areas denote positive anomalies of TEC− M and negative anomalies of M−TEC, respectively.The blue dashed rectangular denotes the storm signals on 26-28 August 2018.F10.7 generally lies between 67-74 sfu (solar flux unit; 1 sfu = 10 −22 W m −2 Hz −1) and peaks 74 fsu on 24 August.The Dst index displays that it is relative quiet, especially 1-10 August, and there is an intense magnetic storm with a maximum depression of − 175 nT (Kp 7 + ) on 26 August.The LTT plots reveal that the EIA (equatorial ionization crest) becomes most intense and moves very poleward on the storm day of 26 August.The most significant positive TEC are observed before those earthquakes on 1,3, 12, 14, 15, 17, and 18 August 2018

Fig. 3
Fig.3Locations of 123 M ≥ 6.0 earthquakes in the Lombok, Indonesia area, where is about ± 15° in latitudes/longitudes from the 5 August 2018 M7.0 earthquake during 2007-2017.The yellow star and red circles denote the M7.0 earthquake (8.3°N, 116.4°E) and the M ≥ 6.0 earthquakes, respectively.The blue open triangle (7.5°N, 115.0°E), near the M7.0 earthquake, is the location of the TEC value extracted from the GIM for the statistical analysis of finding the PEIA characteristics

Fig. 5
Fig. 5 ROC curves for alarming earthquakes based on precursory information indicated in Fig. 4. Color and red curves denote the ROC curves of the observations and the 95% interval, respectively.The red and blue asterisks denote the best point yielding the maximum R score (= TPR−FPR), which is called the Youden index (Youden 1950).p-value is "zero"

Fig. 6
Fig. 6 Distributions of the significant positive TEC anomaly occurrences during the PEIA period associated with the M7.0 earthquake (1-3 August 2018), the M6.s earthquakes (12, 14, 15, 17, and 18 August 2018), and the M7.0 + M6.s earthquakes.For each lattice, the proportion is the number of positive anomalies divided by the total timepoints of the period in percentage.a Overall proportion associated with the M7.0 earthquake.The top proportion associated with b the M7.0 earthquake of 65%, c the M6.s earthquakes of 60%, and d the M7.0 + M6.s earthquakes of 64%.e The overall and f top 75% proportions during the storm period of 26-28 August 2018.The thin black curve denotes the magnetic equator.The heavy black curves present the mid-latitude troughs in the northern and southern hemisphere.Anomalies around ± 65 o N are most likely related to the mid-latitude trough.The black rectangular frame (25°S-5°N, 110°E-150°E) denotes the epicenter area of the M7.0 and M6.s earthquakes

Fig. 7
Fig. 7 Global distributions of significant anomalies at post-noon and post-midnight a-d Ne/Ni, e-h Te/Ti.Red stars denote the epicenter of the 5 August 2018 M7.0 earthquake.The top count of the anomalies associated with i the M7.0 earthquake, j M6.s earthquakes, and k the M7.0 + M6.s earthquakes.It can be seen that in (k), the anomalies in the post-noon Ne, Ni, Te, and Ti, and post-midnight Ni concurrently appear right above the epicenter.The anomalies in Ne/Ni and Te/Ti are defined as greater than the median, 50%, and lower than the first quartile, 25%, respectively

Table 1 M
> 6.0 Earthquakes in Indonesia during August 2018