Seismo-Traveling Ionospheric Disturbances Triggered by the 12 May 2008 M 8 . 0 Wenchuan Earthquake

A network of 6 ground-based GPS receivers in East Asia was employed to study seismo-traveling ionospheric disturbances (STIDs) triggered by an M 8.0 earthquake which occurred at Wenchuan on 12 May 2008. The network detected 5 STIDs on the south side of the epicenter area. A study on the distances of the detected STIDs to the epicenter versus their associated traveling times shows that the horizontal speed is about 600 m s-1. Applying the circle method, we find that the 5 circles intercept at a point right above the epicenter when the horizontal speed of 600 m s-1 is given. Global searches of the ray-tracing and the beam-forming techniques confirm that the STIDs are induced by vertical motions in the Earth’s surface during the Wenchuan Earthquake.


INTRODUCTION
During earthquake occurrences, vertical motion of the Earth's surface creates mechanical disturbances (acoustic gravity waves or AGWs) in the atmosphere, which propagate into the ionosphere and interact with the ionized gas (hereafter, seismo-traveling ionospheric disturbances; STIDs) (Davies 1990). Using ionograms recorded by MF/HF (median/high frequency) ionosondes and frequency shifts probed by HF Doppler sounders, scientists have observed STIDs triggered by strong earthquakes (Davies and Baker 1965;Leonard and Barnes 1965;Row 1967;Yuen et al. 1969;Tanaka et al. 1984;Blanc 1985;Artru et al. 2004;Liu et al. 2006a). Due to limited numbers and/or short distances among stations of ionosonde or Doppler sounding networks, it is rather difficult to study the propagation of STIDs in detail. Liu et al. (2006a) for the first time applied the circle method (Lay and Wallace 1995) on ionograms recorded by 3 ionosonde stations and found the origin and propagation speed of STIDs triggered by the 26 December 2004 M 9.3 Sumatra Earthquake.
Recently, the total electron content (TEC) derived from data recorded by dense ground-based receivers of the global positioning system (GPS) have been employed to examine STIDs such as Rayleigh, gravity, shock, and tsunami waves triggered by earthquakes (Calais and Minster 1995;Afraimovich et al. 2001;Ducic et al. 2003;Artru et al. 2005;Heki and Ping 2005;Astafyeva and Afraimovich 2006;Jung et al. 2006;Liu et al. 2006b;Astafyeva and Heki 2009;Astafyeva et al. 2009;Liu et al. 2010Liu et al. , 2011. Based on Calais and Minster (1995), most scientists find the timedistance relationship between the triggered disturbances and the epicenter to estimate STID speeds. On the other hand, Jung et al. (2006) and Liu et al. (2006b and2010) employed the beam-forming and/or the ray-tracing methods analyzing GPS TEC observations to locate the origin and compute average propagation speeds of STIDs triggered by the 13 January 2001 M w 7.6 El Salvador Earthquake, the 26 December 2004 M w 9.3 Sumatra Earthquake, and the 20 September 1999 M w 7.6 Chi-Chi Earthquake, respectively. The paper briefly reviews the amplification effect, the ground-based GPS TEC observation, and the exiting methods of STID analyses. In data analyses, we first employ the time-distance relationship (Calais and Minster 1995) estimating the propagation speed, and apply the circle method (Lay and Wallace 1995) locating the source, and finally use the ray-tracing (Lee and Stewart 1981) and beam-forming (Huang et al. 1999) method to simultaneously find the speed, source location, and onset time of STIDs triggered by the 12 May 2008 M w 8.0 Wenchuan Earthquake.

AMPLIFICATION FACTOR AND IONOSPHERIC GPS TEC
Near-surface earthquakes causing large vertical nearfield displacement of the Earth's surface excite mechanical disturbances (i.e., AGWs) in the atmosphere, which propagate to the ionosphere where they couple into the ionized gas inducing disturbances by the effects resulting from collision, i.e., STIDs (Blanc 1985). Since atmospheric density decreases almost exponentially with altitude, energy conservation implies that pulse amplitude increases upward as it propagates into the atmosphere and ionosphere. The amplification factor can reach tens to hundreds of thousands at ionospheric heights (Calais and Minster 1995). Both potential and kinetic energy density can be found. Based on the concept of energy conservation, the energy density of STIDs can be expressed as, where G and I denote the ground and the ionosphere, respectively. Thus, G t and v G ( I t and v I ) stand for the mass density of the atmosphere and maximum oscillation speed on the ground (in the ionosphere), respectively. For a simple harmonic motion, the oscillation speed can be written v = Aω, where A is the oscillation amplitude and ω is the angular frequency. Therefore, Eq. (1) can be given as It has been known that when an earthquake occurs, the associated vertical motion of the Earth's surface excites AGWs at numerous frequencies in the atmosphere near the Earth's surface. Since the atmosphere acts a low-pass filter, high frequency waves with periods shorter than 1 minute can quickly be damped. By contrast, low frequency waves with periods longer than 5 minutes generally could travel through the atmosphere efficiently and reach the ionosphere. If the wave-wave interaction is small and/or negligible, the angular frequency of the long-period waves remains unchanged, ω = ω G = ω I , and we find that the oscillation amplitude is inversely proportional to the density, which is written as It is found the density ratio between the ground and the ionosphere in the solar minimum year of 2008 is 10 10 (Kelly 2009). Therefore, for an ideal case of a plane source, the amplification factor is about 10 5 . However, for a real-world case of a point source of an epicenter or a line source of a fault zone, the amplification factor becomes about 10 4 . The GPS consists of more than 24 satellites, distributed over 6 orbits encircling the globe. Recently, geodetic scientists investigated Earth's surface deformation rates using ground-based GPS networks (e.g., see Calais and Amarjargal 2000). While observing deformation, the same networks can be simultaneously employed to monitor the ionospheric TEC (see papers listed therein Liu et al. 1996Liu et al. , 2001. To simultaneously monitor a large area of the ionosphere, a network of ground-based receivers GPS is ideal to be used. Each satellite transmits two frequency signals, 1575.42 and 1227.60 MHz. Since the ionosphere is a dispersive medium, scientists are able to evaluate the ionospheric effects on the radiowave propagation or the corresponding ray path TEC with measurements of the modulations on carrier phases and code pseudoranges recorded by dual-frequency receivers (cf., Liu et al. 1996). The TEC is defined as the sum (or integration) of the ionospheric electron density along a line-of-sight from a ground-based receiver to the associated GPS satellite at an orbital altitude of about 20200 km (for detail, see Liu et al. 1996). Here, the line-of-sight or slant TEC (STEC) between a GPS satellite and a ground-based receiver can be written as where L 1 and L 2 in meter denote the carrier phases of the two frequencies f 1 and f 2 in hertz, and δr and δs are the differential biases for receiver and satellite, respectively. The greatest electron density in the ionosphere usually situates at about a 300-km altitude, which contributes to the heaviest weight in the slant TEC calculation. Based on the spirit of the center of mass, the center of the TEC at about a 350-km altitude is termed the ionospheric point. Therefore, a GPS TEC acts as a space seismometer floating at about 350 km above the Earth's surface and monitors ionospheric disturbances. From recorded GPS broadcast ephemeris, the latitudinal and longitudinal coordinates of each space seismometer can be computed. Assume the ionosphere to be a thin shell at 350 km altitude, the slant TEC along the ray path can be converted, usually using a slant function of the satellite zenith, into the vertical TEC (VTEC, for simplicity hereafter, TEC) at the space seismometer location (Tsai and Liu 1999).

SEARCH THE STID ORIGIN
The method of intersecting circles might be the first employed for locating the hypocenter (source, or epicenter) of an earthquake (Lay and Wallace 1995). Scientists calculate the possible circular distance (zone) to the source of seismic waves from each station, which is equal to the product of the propagation speed and the traveling time. When circular zones of three or more stations intersect (i.e., triangulation) at the same location, we then consider the source to be located. The short coming of this method is that the onset time of the earthquake and the propagation speed of the seismic waves should be roughly known in advance.
On the other hand, the ray-tracing technique (Lee and Stewart 1981) and the beam-forming technique (Huang et al. 1999) also have been employed traditionally by seismologists to locate an earthquake source (or hypocenter). In the ray-tracing technique, for a given velocity model, scientists estimate the location of a hypocenter and calculate the arrival time of seismic waves at each seismometer (the forward calculation). From the difference between the calculated and observed arrival times at the seismometers, a new location of the hypocenter is further issued (or given) and tried. When the differences reach the minimum, the iteration process is stopped, and the estimated center is then declared to be the hypocenter. Similarly, scientists can guess the location of a hypocenter but calculate the earthquake onset time from each seismometer (using a backward calculation). The minimum in differences of the calculated onset times is the criterion of the hypocenter location determination. Here, we consider the STID as a seismic wave recorded by a space seismometer orbiting at an altitude of 350 km, and adopt the backward calculation. Assume the average speeds in the vertical and horizontal direction to be V Z and V H , respectively. The traveling time to each STID point (i.e., space seismometer "i"), which is the sum of traveling times in the vertical and horizontal directions, can be written as, where the STID at the ionospheric point altitude ΔZ i, j = 350 km and the horizontal distance from the STID to the trial epicenter ΔS i, j . Since the arrival time t i at a STID is recorded, the onset time t Gi at the trial epicenter "j" can be computed, If N points of STIDs are detected, the standard deviation of the computed onset times at the trial center can be written as where μ j is the associated mean value. The minimum in the standard deviation is the criterion of the location determination of the STID origin/source. The average onset time of the minimum in the standard deviation is used to find when and where the STIDs are activated. By contrast, for a given onset time, the beam-forming technique (Huang et al. 1999) estimates a hypocenter and computes the speeds by dividing the distances from the hypocenter to the seismometers by the differences between the onset time and observed arrival times at the seismometers. Similarly, when the standard deviation of the computed speeds yields the minimum, the hypocenter is considered to be located. Here we simply guess a onset time t o and compute the velocity V i, j by dividing the distance from the STID point to the trial epicenter ΔS i, j by a given traveling time Δt i (= t i -t o ) at each STID point, Similarly, for N points of the detected STIDs, the standard deviation of the computed velocities (speeds) at the trial center can be written as The minimum in the standard deviation is the criterion of the location determination of the STID origin/source. The average speed of the minimum in the standard deviation is considered the STID propagation.

OBSERVATION AND INTERPRETATION
A magnitude M 8.0 earthquake occurred in Wenchuan, China (31.0°N, 103.4°E geographic; 24.9°N, 175.5°E geomagnetic) with a depth of 19 km at 0628 UT (universal time) on 12 May 2008. Figure 1 displays the Wenchuan epicenter as well as the locations of 6 GPS receivers and their associated fields of coverage/observation. A high-pass filter with cutoff period less than 10 minutes is applied to analyze the TECs observed at the 6 stations. It can be seen that 5 pronounced fluctuations of STIDs observed at the GPS station kunm at 0635 -0720 UT which was about 10 minutes after the earthquake occurrence. Figure 2 illustrates temporal variations of the filtered TEC and the associated distance to the epicenter versus time. A line of the least square fitted of the 5 STID arrival times intersects at 0636 UT. Since the reported earthquake onset time is 0628 UT, it takes 480 sec (= 8 min) for the STID traveling into the ionosphere. Therefore, the average horizontal and vertical speeds are about 557.2 ± 50 m s -1 and 729 (= 350 km/480 sec) m s -1 , respectively.
Since the STIDs are assumed to be detected at a fixed altitude of 350 km, we can simply consider their horizontal distances and speeds. Thus, for the circular method, there-fore we let the onset time be 0636 UT instead of 0628 UT. Various horizontal speeds 400 -800 m s -1 have been tested. Figure 3 illustrates that when the horizontal speed of 600 m s -1 is given, the 5 circles intersect at one location, where is right above (or slightly west of) the Wenchuan epicenter.
In both the ray-tracing and beam-forming calculations, we have tested the trial center within 20 -33°N and 95 - Fig. 1. Locations of the Wenchuan epicenter and 6 ground-based GPD receivers. TEC filtered with a high-pass filter with a cutoff period less than 10 minutes. Pronounced fluctuations of 5 STID signatures were observed by the GPS station kunm. The earthquake magnitude M 7.9 reported by US glological survery (USGS) is later ranked as M 8.0 by China Earthquake Administration (CSA). The red star denotes the epicenter. 111°E shifting by 0.1°. For the ray-tracing search, we try with the horizontal speed V H ranging from 200 to 1000 m s -1 . The optimal results can be obtained by finding the minimum values of the standard deviation of the calculated onset times, the time difference between the calculated average and reported onset times, as well as the distance between the calculated center (disturbance origin) and the Wenchuan epicenter. It is found that when the horizontal speed 620 m s -1 , the contour converges at 30.4°N, 102.9°E where is about 86-km southwest of the Wenchuan epicenter (Fig. 4). Note that the Wenchuan epicenter is within the standard deviation of 25 seconds. The associated onset time is 0628 UT (= 0635 UT-7 minutes of the vertical traveling time).
For the beam-forming technique, we try with the onset time varying from 0600 to 0800 UT and compute the mean speed and standard deviation from each trial center to the 5 STIDs. Figure 5 displays that the onset time at 0638 UT yields the best result (the minimum of the standard deviation) of the STIDs origin being located at 31. 1°N, 103.0°E where is about 30 km west of the Wenchuan epicenter. The optimal result of the mean horizontal speed and standard deviation is V H = 614 ± 15 m s -1 traveling away from the STID origin. Liu et al. (2006a) observed 2 types of STIDs with traveling speeds of 2.5 km s -1 and 370 m s -1 triggered by the 24 December 2004 M 9.3 Sumatra earthquake. Liu et al. (2010) found that the STIDs triggered by the M 7.6 Chi-Chi earthquake range from 700 to 900 m s -1 . Astafyeva et al. (2009) analyzed STIDs induced by the 4 October 1994 M 8.1 Kurile earthquake and found that starting from about 600 -700 km away from the epicenter, the disturbance seems to divide into two perturbations which propagate with different ve-locities, one at about 3 km s -1 and the other at about 600 m s -1 and suggest that the former perturbation is triggered by acoustic waves launched by a Rayleigh surface wave and the later by the vertical displacements of the ground motion due to the earthquake. Figure 2 shows that the distances of the 5 STIDs to the epicenter are greater than 540 km (= 600 m s -1 × 15 min). Artru et al. (2004) and Heki and Ping (2005) show that the acoustic speed is about 300 m s -1 from the Earth's surface to the mesosphere at about 90 -100 km altitude, and speed up to 1000 -1100 m s -1 at 300 km altitude in the ionosphere. Therefore, if the triggered disturbances depart with low elevation angles, they travel nearly horizontally between the troposphere and the mesosphere with speeds of about 300 m s -1 , and then both locally and vertically perturb the ionospheric TEC (Liu et al. 2006a). Notably, those disturbances depart with high elevation angles, travel nearly vertically into the thermosphere (ionosphere), and propagate horizontally interacting with the ionized gas, which should travel with average speeds of 600 -1100 m s -1 (Astafyeva et al. 2009;Liu et al. 2010). Based on the above observations and simulations, we speculate that large vertical motions of the Earth's surface near the Wenchuan epicenter excited AGWs in the atmosphere, which propagated obliquely with a moderate elevation angle into the ionosphere as STIDs traveling at 600 m s -1 . Although the distance from the epicenter to the GPS station wuhn is not so different, the 5 STIDs are observed only by the GPS station kunm. Scientists hypothesize that the northward (Heki and Ping 2005;Otsuka et al. 2006), and eastward/westward propagating disturbances might be attenuated selectively by interaction between the movements of charged particles in STIDs and the Earth's magnetic fields. Note that the STIDs propagations with a moderate elevation angle in the northward, eastward, and westward directions are perpendicular to the magnetic fields. This once again suggests that AGWs near the Wenchuan epicenter obliquely propagate with moderate elevation angles into the ionosphere. It is very difficult to obtain ground-based GPS measurements from China. Because no data is available near the epicenter, STIDs induced by the Rayleigh waves cannot be found or identified. Although the data are very limited, the timedistance relationship, the circle method, and the ray-tracing and beam-forming techniques confirm that the STIDs travel with an average speed of 600 m s -1 and were induced by the Wenchuan earthquake.   (a) (b) (c)