An observing system simulation experiment for FORMOSAT-5/AIP detecting seismo-ionospheric precursors

AIP (advanced ionospheric probe) on board FORMOSAT-5 is a next version and advanced generation science payload of IPEI (ionospheric Plasma and Electrodynamics Instrument) on board ROCSAT-1 (i.e., FORMOSAT-1). We examine the ion density, ion temperature, and ion velocity probed by ROCSAT-1/IPEI, as well as the global ionospheric map (GIM) of the total electron content (TEC) derived by ground-based GPS receivers during the 31 March 2002 M6.8 Earthquake in Taiwan to determine whether FORMOSAT-5/AIP can be used to detect seismo-ionospheric precursors (SIPs). It is found that in the epicenter area, the GIM TEC significantly decreases 1 5 days before the earthquake, while the ROCSAT-1/IPEI ion density significantly decreases and ion velocity in the downward direction anomalously increases. The anomalous decreases in the ROCSAT/IEPI ion density and the GIM TEC concurrently appear around the epicenter area 1 5 days before the earthquake, which suggests that FORMOSAT-5/AIP can be employed to detect SIPs. Moreover, the increase in the downward velocity implies that a westward electric field generated during the preparation earthquake period is essential. Article history: Received 11 April 2016 Revised 27 June 2016 Accepted 18 July 2016

The short comings of these studies are that the observations were taken in very limited areas and mainly over land.Ground-based satellite receivers have recently provided an opportunity for conducting pseudo global observations.Liu et al. (2001) pioneered employing the total electron content (TEC) derived from a local network of ground-based GPS receivers to detect SIPs appearing 1, 3, and 4 days before the 21 September 1999 M7.6 Chi-Chi earthquake.Similarly, Liu et al. (2004a) statistically studied the GPS TEC and 20M ≥ 6.0 earthquakes in Taiwan from September 1999 to December 2002.They also reported that the GPS TEC significantly decreased 1 -5 before the earthquakes.Liu et al. (2004b) once again statistically investigated variations in the foF2/ TEC and 144M ≥ 5.0 earthquakes in Taiwan during 1997 -1999.They confirmed that the foF2 and TEC yielded similar tendencies, and often concurrently registered pronounced decrease anomalies 4 days before the earthquakes.Following these discoveries, the global ionosphere map (GIM) of the GPS TEC has been used to globally study SIPs appearing a few days before the 26 December 2004 M9.3 Sumatra-Andaman Earthquake, the 11 May 2008 M8.0 Wenchuan earthquake and the 12 January 2010 M7.0 Haiti earthquake (Liu et al. 2009(Liu et al. , 2010a(Liu et al. , 2011)).Nevertheless, the GIM TEC observations still suffer from limited ground-based receivers over the ocean and remote areas.
Satellites generally provide global, uniform coverage Terr. Atmos. Ocean. Sci., Vol. 28, No. 2, 117-127, April 2017 for studying SIPs (Pulinets 1998;Afonin et al. 1999;Hayakawa et al. 2000).DEMETER (Detection of Electro-Magnetic Emissions Transmitted from Earthquake Regions) might be the first satellite designed specifically to detect SIPs.It is a micro-satellite that observes the region within 65° geomagnetic latitude at a 680 km altitude, 98.3° inclination orbit during June 2004 to December 2010.Its main scientific objectives are the detection and characterization of ionospheric electrical and magnetic disturbances in connection with seismic activity (Cussac et al. 2006).DEMETER carries 6 payloads that detect the plasma density and ion composition (Berthelier et al. 2006); measuring the electron density and temperature (Lebreton et al. 2006); performing an automatic identification and classification of the whistlers from ELF -VLF electric field measurements (Elie et al. 1999), and records the electric field, magnetic field and high energy particles, respectively.Many research articles reported density variations and amplitude changes in electro-magnetic emissions above earthquake regions using DEMETER observations (cf.http://smsc.cnes.fr/DEMETER/A_ publications.htm).
Although many ionospheric scalar quantities, especially, the TEC, electron density, electron temperature, ion density, ion temperature, etc., for SIPs have been investigated intensively, the vector quantity, such as the ion velocity, has not yet been reported in detail.ROCSAT-1 is a low-earthorbit scientific experimental satellite with a circular orbit and a 35-degree inclination at 600 km altitude during January 27 1999 to 17 June 2004.Its main payload, IPEI (ionospheric Plasma and Electrodynamics Instrument) provides the most comprehensive ion density, the ion temperature, and ionospheric velocity observations at the satellite orbit (Yeh et al. 1999).Advanced Ionospheric Probe (AIP), a next generation of IPEI, is a FORMOSAT-5 science payload, to be launched in a sun synchronous orbit at 720-km altitude with 98.28-degree inclination angle, 99.19-minute period, a local time of descending (ascending) node at 9: 45 -10:15 (21:45 -22-15) LT (local time), a revisit time of every other day.This satellite will be launched in 2016.AIP is an all-in-one plasma sensor with a sampling rate up to 8192 Hz to measure the ion density, ion velocities, ion temperature, electron temperature, and ambient magnetic fields at the satellite orbit over a wide range of spatial scales.With the comprehensive data set available from AIP, scientists are able to conduct a systematic examination of the longitudinal, latitudinal, and seasonal variations in the topside F-region plasma parameters.The FORMOSAT-5's orbit allows us to construct daily daytime and nighttime global maps of each AIP plasma parameter with a longitudinal resolution of 12-degrees at global fixed local times of about 10:00 LT and 22:00 LT, respectively.Although simply mapping at the two global fixed local times, FORMOSAT/AIP yields an advantage of observing EIA (equatorial ionization anomaly) crests (cf.Lin et al. 2007), middle latitude trough (cf. Lee et al. 2011), etc. in the Northern and Southern hemispheres in the same orbit.The data are designed for monitoring the ionospheric space weather and also for detecting precursors (i.e., SIPs) and co-seismic disturbances (i.e., seismotraveling ionospheric disturbances) induced by worldwide earthquakes (cf.Liu et al. 2016).
This paper presents an observing system simulation experiment of ROCSAT-1/IPEI searching SIPs associated with the M6.8 earthquake occurring at (24.3°N, 122.2°E) and 14:52:40 LT on 31 March 2002 (hereafter the 331 earthquake).This experiment is used to determine if FOR-MOSAT-5/AIP can be used to detect SIPs associated with earthquakes.While, the GIM (ftp://cddisa.gsfc.nasa.gov/pub/gps/products/ionex) published at the Center for Orbit Determination in Europe (CODE) of the TEC searches the global distribution of SIPs for the 331 earthquake, the ion density, ion temperature, and ion velocity concurrently measured by IPEI are employed to detect SIPs.We report and discuss for the first time how ion velocity SIPs could lead to better understanding of possible earthquake causal mechanisms.

ObSeRvATIOnS And dATA AnAlySIS
The CODE GIM of the TEC maps, with 2-h time resolution (12 time points in a day) covering ±87.5°N latitude and ±180°E longitude ranges with spatial resolutions of 2.5° and 5°, respectively, is ideal for the detection of temporal and spatial SIPs associated with the 331 earthquake.Each map consists of 5183 (= 71 × 73) grid points (lattices).A quartile-based process is performed to detect GPS TEC variation anomalies.At each time point, we compute the median M of every successive GPS TEC 15 days as well as the deviation between the observation on the 16 th day and the computed median.We also calculate the first (or lower) and the third (or upper) quartiles, denoted by LQ and UQ, respectively to provide the deviation information.Note that assuming a normal distribution with mean m and standard deviation σ for the GPS TEC, the expected LQ and UQ values should be m -0.67σ and m + 0.67σ, respectively (Kotz et al. 2005).We set the lower bound, LB = M -1.5(M -LQ) and upper bound, UB = M + 1.5(UQ -M) to have stringent criterion.Therefore, the probability of a new GPS TEC in the interval (LB, UB) is approximately 68%.The median together with the associated LB and UB then provides references for the GPS TEC variations on the 16 th day.Thus, when an observed GPS TEC on the 16 th day is not in the associated (LB, UB), we declare an upper or lower anomalous GPS TEC signal.Figure 1 shows a GIM TEC time series together with its associated median and upper/lower bound above the epicenter isolated from the GIM maps during 16 March to 14 April 2002.Since the statistical results show that the GPS TEC significantly decreases 1 -5 before large earthquakes in Taiwan (Liu et al. 2000(Liu et al. , 2004a(Liu et al. , b, 2006)), the anomalous TEC decreases, exceeding the associated lower bound, appearing on 1, 2, 4, and 5 days before, which are possibly related to the 331 earthquake.The longest period appears from 20:00 UT on 28 March to 14:00 UT on 29 March 2002.
In addition to studying temporal TEC variations over the epicenter, we also examined TEC variations along the Taiwan longitude 121°E.Figure 2 illustrates the latitudetime-TEC (LTT) plots along the Taiwan longitude extracted from the GIM between 15 days before and 1 day after the earthquake (16 March to 1 April).It can be seen that during the study period the two EIA crests become weak and tend to move equator-ward 2 -7 days before the 331 earthquake.The two EIA crests reaching the minimum values move sig-nificantly equator-ward on 29 March, in coincidence with the TEC over the epicenter decreasing extremely on the same day (Fig. 1).These EIA crest anomalous features imply that the daytime eastward electric field of the daily dynamo is being significantly disturbed and reduced (Kelley 2009).
The GIM TEC over Taiwan and along the Taiwan longitude decreases significantly and yields the minimum value on 29 March 2002.To discriminate the local effect of earthquakes from global effects such as solar disturbances, magnetic storms, etc. and to see if the TEC decrease anomalies appear specifically in the earthquake region for a long duration or simply randomly occur worldwide on 29 March 2002, a spatial analysis on the global distribution of decrease anomalies in the GIM TEC is required.Figure 3 The decrease anomaly lattices with 1 (00:00 -00:00 UT), 2 (00:00 -02:00 UT), 3 (00:00 -00:04 UT), and more repeating time points (times) occur in low-latitudes worldwide, whereas those with 4, 5, 6, and more repeating times appear mainly in the epicenter area, and its magnetic conjugate area and the South American area.It can be seen the decrease anomaly lattices repeat up to 7, 8, 9, and even more time points that appear repeatedly mainly around the epicenter area.Those around the conjugate point and the South American area disappear drastically.The decrease anomaly lattices with 10 -11 and more repeats intersect/appear in two small areas around the South Atlantic Anomaly and south of the epicenter.Finally, the decrease anomalies appearing specifically in a small area slightly south of the epicenter, repeat up to 12 times (i.e., 24 = 2 × 12 hour duration) on 29 March 2002 (remember, the GIM with 2-h time resolution; 12 time points of a day).Therefore, the decrease anoma-lies in the GIM TEC appear persistently in the small area slightly south of the epicenter, for the entire universal day of 29 March 2002.The long-persisting anomalies on 29 March indicate that these anomalies are most likely related to the 331 earthquake.
Since Figs. 1 and 2 reveal that the GIM TECs significantly decrease in the afternoon period of 26 -30 March 2002, we then further examined the global distribution of occurrences (or counts) of the decrease anomalies appearing in the afternoon period at 4 time points, 12:00 -18:00 LT (06:00 -12:00 UT), during 1 -5 days before the 331 earthquake (Fig. 4).In total, there are 20 (= 4 × 5) time points in the 5-day period.A lattice at 20°N, 115°E, which is again located slightly south of the epicenter, yields the maximum occurrence of 17 counts (85% = 17/20).Figures 3 and 4 show that the most frequent decrease anomalies fall in an area (18 -23°N, 110 -130°E) denoted by a black rectangular frame.
Figures 5a, b, and c illustrate the ion density observed in the afternoon on Taiwan (i.e., 06:00 -12:00 UT) 1 -5 Fig. 3. Locations of the decrease anomaly in the GIM TEC repeatedly appear at various time points on the universal day of 29 March 2002.The number or duration of the repeat time point is noted at the top of each panel.The red lattices denote the decrease anomaly intersection.The star symbols denote the epicenter.The UT 00 -00, 00 -02, 00 -04, etc. stand for the negative anomalies continuously appear over the time period of 00:00 -00:00, 00:00 -02:00, 00:00 -04:00 UT, etc., respectively.(Color online only) days before the 331 earthquake, 45 days before/after (as a reference) the 331 earthquake, and the difference in corresponding reference subtracted from the 5-day observations, respectively.Figure 5c displays that the ion density around the epicenter area denoted by the black rectangular frame (18 -23°N, 110 -130°E) significantly decreases in the 5-day period, which agrees with the GIM TEC decreases shown in Figs. 3 and 4. We further statistically investigate the ion density anomalies with box-and-whisker (box) plots (Wilcox 2010) within the rectangular area in the 5-day (Fig. 5a) and in the reference (Fig. 5b).The lower (upper) quartile is the number such that at least 25% (75%) of the sorted observations are less (greater) than or equal to it.In the plot the ends of the box are the upper and lower quartiles.The horizontal line in the box denotes the statistical median.The dots are the outliers that exceed 1.5 times the inter-quartile range from the end of the box and are unusually small or large.The lines extending out from the box are called whiskers.The ends of the whiskers mark the smallest and largest values that are not outliers.When one median (i.e., horizontal line) does not overlap with the other's upper and lower quartiles (i.e., box), the two are declared to be significantly different.Figure 5d displays that the two boxes do not overlap each other at all, which can be declared different under significance level 0.06 (= 0.25 × 0.25) (Wilcox 2010).This confirms that the ion density over the epicenter area significantly decreased 1 -5 days before the 331 earthquake.
Similar analyses were also applied to the ion temperature and the ion velocity.Figures 6a, b, and c depict within the rectangular, no clear difference in the ion temperature between the 5-day and the associated reference.Figure 6d confirms that the two boxes overlap each other, which indicates that no SIPs of the ion temperature related to the 331 earthquake can be detected.Figures 7a -i illustrate the 5-day observation, the reference, the difference in the ion velocity parallel to the Earth's magnetic field (V par ); perpendicular to the Earth's magnetic field in the radical upward component (V PerM , positive in the upward direction); and perpendicular to the Earth's magnetic field in the zonal component (V PerZ , positive in the eastward direction), respectively.Figure 7j shows that V par and V PerM significantly decrease 1 -5 days before the 331 earthquake, while V PerZ reveals no obvious difference.It can be seen that the median values in V par (V PerM ) during the 5-day and the reference days are -21 and 40 m s -1 (-10 and 18 m s -1 ), respectively.Thus, the ion velocity in the anti-parallel magnetic field direction (i.e., equator-ward) and the perpendicular magnetic field downward direction significantly increase.Kamogawa et al. (2004) reported that the TEC derived from a local ground-based GPS receiver network significantly decreased 2, 4, and 5 days before the 331 earthquake.A corona probe near the epicenter recorded that the atmospheric electric field anomalously fluctuates 1 day and 40+ minutes before the earthquake.Statistical results show that the ionospheric electron density and the GPS TEC tend to decrease significantly 1 -5 days before large earthquakes in Taiwan (Liu et al. 2000(Liu et al. , 2001(Liu et al. , 2004a(Liu et al. , b, 2006)).Figures 1 and 2 reveal that the TEC anomalously decreases over the epicenter and the EIA crests 1, 2, 4, and 5 days before the 331 earthquake, which agrees well with the TEC derived from the local GPS network (Kamogawa et al. 2004) and the SIP statistical result in the TEC derived from the local GPS network in Taiwan (Liu et al. 2004a).

dIScuSSIOn And cOncluSIOn
Figures 3 and 4 display that the GIM TEC decrease anomalies appear specifically in the small area south of the epicenter 1 -5 days before, especially on 29 March 2002 (2 days before), the earthquake.Figure 5c also illustrates that the ROCSAT-1/IPEI ion density specifically decrease in the small area south of the epicenter 1 -5 days before the 331 earthquake.The GIM TEC decrease anomalies and the ROCSAT-1/IPEI ion density were specifically and concurrently detected in the small area south of the epicenter 1 -5 days before, indicating that the SIPs are related to the 331 earthquake.
Although Fig. 6 shows no clear SIPs in the ROC-SAT-1/IPEI ion temperature before the 331 earthquake, Ho et al. (2013a, b) and Liu et al. (2015) found that the DEMETER electron density, ion density and ion temperature were useful/sensitive for detecting the 2010 M8.8 Chile earthquake and the 2008 M8.0 Wenchuan earthquake, while Oyama et al. (2011) reported possible precursor features in the ionospheric ion density associated with the 16 October 1981 M7.5 earthquake (-33.13°N, 73.07°E).Therefore, the ionospheric ion temperature would still be a physical quantity for detecting SIPs.Despite no electron temperature available from ROCSAT-1/IPEI for the observing system simulation experiment of FORMOSAT-5/AIP, we find that Oyama et al. (2008) examined ionospheric electron temperatures observed by the HINOTORI satellite during three earthquakes; M6.6 occurred in November 1981, M7.4 and M6.6 in January 1982 over the Philippines.They clearly show that the electron temperature around the epicenters significantly decreased in the afternoon periods within 5 days before and after the three earthquakes.Thus, the FOR-MOSAT-5/AIP electron temperature would be also useful in observing SIPs in the future.Liu et al. (2010b) statistically studied the ionospheric EIA crests of the GPS TEC response to all 150M ≥ 5.0 earthquakes in the Taiwan area during 2001 -2007.They found that along the Taiwan longitude the EIA crest significantly moves equator-ward and appears at an earlier time in the afternoon period 1 -5 days before the earthquakes.These results suggest that seismo-atmospheric electric (westward) fields can be generated around the epicenter of a forthcoming earthquake, which map along the Earth's magnetic field line to the low-latitude and/or equatorial ionosphere, and then weaken the daytime daily dynamo fountain (eastward) electric field during the preparation period.Figure 2 depicts that the two EIA crests become weak and tend to move equator-ward 2 -5 days before the 331 earthquake, which agrees well with Liu et al. (2010b).This agreement suggests that seismo-generated westward electric fields appear 1 -5 days before the 331 earthquake.
The EIA is a typical F-region ionosphere phenomenon at low latitudes.It is characterized by an electron density trough at the magnetic dip equator and a dual band of enhanced electron density (two crests) at about 15° north and south of the trough.The EIA formation is a result of the diurnal variation in the zonal electric field, which primarily points eastward during the day.With the horizontal geomagnetic field at equatorial latitudes, the plasma is lifted up by a vertical e × b drift.Once transported to higher altitudes, the plasma diffuses downward along the geomagnetic field lines into both hemispheres due to gravitational and pressure gradient forces.The overall process is known as the ''e × b plasma fountain effect'' (for detail see, Stolle et al. 2008).Thus, a weaker fountain should result in a more equator-ward motion and a less dense of the crest.
Figure 7 reveals that the speeds in V par (V PerM ) around the epicenter during the 5-day and the reference days are about -21.43 and 40.41 m s -1 (-11.12 and 17.05 m s -1 ), respectively.Based on the dynamo theory, the electric field can be expressed as e = -v × b, where v is the ion velocity and the b is the Earth's magnetic field.The seismo-generated electric field during the SIP days can be estimated e sg = e SIP -e 0 , where e SIP and e 0 stand for the electric fields on the precursor day and the reference day, respectively.From the IGRF (international geomagnetic reference field) model, we find the b field at 600 km altitude over the epicenter is 0.32 × 10 -4 T with the magnetic dip I = 33.4degree.No electric field would be produced since V par is parallel to the magnetic field.By contrast, for the SIP day, V PerM = -11.12m s -1 , e SIP = 0.36 mV m -1 in the westward direction while for the reference day, V PerM = 17.05 m s -1 , e 0 = 0.55 mV m -1 in the eastward direction.Therefore, the seismo-generated electric field can be calculated using e sg = e SIP -e 0 = 0.91 mV m -1 in the westward direction.Based on Mozer and Serlin (1969) and Kelley (2009), the atmospheric field can be non-attenuatedly mapped along the same magnetic field in the atmosphere, ionosphere and magnetosphere, and thus the atmospheric electric field in the magnetic westward direction near the Earth's surface should be about e sg_EW = 0.91 mV m -1 .
Because V par is parallel to the magnetic field, no electric field would be produced, and we might be able to roughly estimate the parallel seismo-generated electric field e sg_par = e SIP_par -e 0_par using electric field mapping from the near Earth's surface to the ionosphere (Mozer and Serlin 1969;Kelley 2009).Based on Kelley (2009), the field aligned current J par = J i_par + J e_par = nqv par + J e_par = σ 0 e par , where J i_par and J e_par stand for the field aligned ion and electron current, while n, q, σ 0 , and e par are the ion (or electron) density, the charge of 1.6 × 10 -19 C, the direct (i.e., aligned) conductivity, and the field aligned electric field, respectively.The direct conductive σ 0 = 60 mho m -1 (Kelley 2009), while the results in this study reveal that the ion density and the ion velocity during the 5-day SIP days (the reference days) are 9.7 × 10 11 # m -3 and -21.43 m s -1 (1.35 × 10 12 # m -3 and 40.41 m s -1 ), respectively.Since there is no data for J e_par (or the electro drift velocity), the parallel seismo-generated electric field is simply approximated using e sg_par = J i_SIP_par / σ 0 -J i_0_par / σ 0 = (nqv par ) _SIP / σ 0 -(nqv par ) _0 / σ 0 = [9.7 × 10 11 × 1.6 × 10 -19 × (-21.43)]/ 60 -(1.35 × 10 12 × 1.6 × 10 -19 × 40.41) / 60 = -5.43 × 10 -8 -1.44 × 10 -7 = -1.5 × 10 -7 (V m -1 ) = -0.15(μV m -1 ).That is the parallel seismo-generated electric field is about 0.15 μV m -1 in the anti-magnetic field direction in the ionosphere.Moreover, if we assume the electric field mapping without attenuation from the ground to the ionosphere, the seismo-generated electric field in the near Earth's atmosphere would be e sg_ver = e sg_par / cos(33.4)= 0.18 μV m -1 in the upward direction.Similarly, we could again estimate the meridian horizontal component of the atmospheric electric field near the Earth's surface as e sg_NS = e sg_ver × sin(33.4)= 0.10 μV m -1 in the magnetic southward direction.Due to the lack of electron velocity measurements, the magnitude could not be obtained correctly, however, the polarity upward and southward of the seismo-generated electric fields were found.Nevertheless, the seismo-generated electric field in the westward direction should be much greater than that in the vertical or the magnetic southward direction.
In conclusion, the ROCSAT-1/IPEI ion density and GIM TEC anomalies appeared specifically over the epicenter 1 -5 days before the 331 earthquake.This agreement together with the previous satellite results (such as, Oyama et al. 2008;Liu et al. 2015) suggests that FORMOSAT-5/ AIP could be employed to detect SIPs of the ion density, ion velocities, ion temperature, and electron temperature.Both GIM TEC and ROCSAT-1/IPEI ion velocities show that the seismo-generated electric field is essential in causing SIPs.This is the first time that the seismo-generated electric field was computed using observation data, which in turn suggests that the ion velocity offers a better understanding of SIP causal mechanisms.It is worth mentioning that although in the case study on the 331 earthquake, a negative effect in the plasma (i.e., ion/electron) density was observed and both negative and positive effects before earthquakes were reported.It would be important to find the positive/negative effects of the plasma density response to the plasma velocity (i.e., the seismo-generated electric field) in the near future.

Fig. 1 .
Fig. 1.The TEC over the epicenter area extracted from GODE GIM during 16 March to 14 April 2002.The red, black dotted, and two black solid curves denote the GIM TEC, associated median and upper/lower bound (UB/LB), respectively.The LB and UB were constructed using the 1 -15 previous days' moving median (M), lower quartile (LQ), and upper quartile (UQ).Here, LB = M -1.5(M -LQ) and UB = M + 1.5(UQ -M).Red and black shaded areas denote differences in O-UB and LB-O, respectively, where O is the observed GIM TEC.(Color online only)