Dependence of waveform of near-field coseismic ionospheric disturbances on focal mechanisms

Using Total Electron Content (TEC) measurements with Global Positioning System we studied ionospheric responses to three large earthquakes that occurred in the Kuril Arc on 04 October 1994, 15 November 2006, and 13 January 2007. These earthquakes have different focal mechanisms, i.e. high-angle reverse, low-angle reverse, and normal faulting, respectively. TEC responses to the 2006 and 2007 events initiated with positive and negative changes, respectively. On the other hand, the initial TEC changes in the 1994 earthquake showed both positive and negative polarities depending on the azimuth around the focal area. Such a variety may reflect differences in coseismic vertical crustal displacements, which are dominated by uplift and subsidence in the 2006 and 2007 events, respectively, but included both in the 1994 event.


Introduction
A large earthquake excites surface waves that go round the Earth several times, and body (P and S) waves that penetrate the core and reach the opposite side of the Earth. Nearly eighty years ago Shida (1929) discovered that polarity of initial motions of P-waves has a clear systematic distribution. In modern seismology, such initial motion information is widely utilized to study focal mechanisms of earthquakes, and helps us study tectonic backgrounds of earthquake occurrences.
Earthquakes also excite atmospheric waves that propagate in the Earth' atmosphere (Blanc, 1985;Artru et al., 2004). Vertical displacements of the ground induce pressure waves in the neutral atmosphere. They propagate upward and grow in amplitude by several orders of magnitude as they attain ionosphere heights through the decreasing atmospheric density. Such waves can initiate the ionospheric plasma motion due to the collision interaction between neutral and charged particles, and produce perceptible perturbations in the ionosphere electron density (Calais and Minster, 1998;Artru et al., 2001;Heki and Ping, 2005;Astafyeva and Afraimovich, 2006;Dautermann et al., 2008). Do the initial motion polarities of such ionospheric disturbances depend, like P-waves, on the earthquake focal mechanisms?
Ionospheric disturbances after large earthquakes (Coseismic Ionospheric Disturbances, CID) are the superposition of perturbations by different origins such as coseismic vertical crustal movements (Calais and Minster, 1995;Afraimovich et al., 2001Afraimovich et al., , 2006Heki and Ping, 2005;Heki et al., 2006), propagating Rayleigh surface waves (Ducic Copyright c The Society of Geomagnetism and Earth, Planetary and Space Sciences (SGEPSS); The Seismological Society of Japan; The Volcanological Society of Japan; The Geodetic Society of Japan; The Japanese Society for Planetary Sciences; TERRAPUB. et al., 2003;Artru et al., 2004). Apart from that, large submarine earthquakes often generate large tsunamis, and they also can produce lower-frequency perturbations in the ionosphere through internal gravity waves (Artru et al., 2005;Liu et al., 2006;Occhipinti et al., 2008). Attributes of these perturbations (e.g., propagation velocities, waveforms, periods, etc.) depend on their origin as well as the distance from the epicenter. For example, ionospheric disturbances excited by Rayleigh waves propagate with a velocity ∼3.3 km/s and become more pronounced as they travel farther due to their small geometric decay. Perturbations induced by propagating tsunamis have quasi-periodic waveforms and travel with a velocity about 190-220 m/s over thousands of kilometers.
Here we focus on near-field CID that is excited by coseismic vertical crustal movements and propagates through the atmosphere as an acoustic wave. Such CID is usually detected in the vicinity of the epicenter (within ∼500 km) 10-15 min after the main shock. They propagate with a velocity ∼800-1000 m/s, faster than the gravity waves but slower than the Rayleigh waves. Their waveforms have been described as "N -type" waves, consisting of leading and trailing shocks connected by smooth linear transition regions (Pavlov, 1986). The parameters of such CID were well examined by Total Electron Content (TEC) measurements with Global Positioning System (GPS) (Calais and Minster, 1995;Afraimovich et al., 2001Afraimovich et al., , 2006Heki and Ping, 2005;Heki et al., 2006;Lognonne et al., 2006).
So far little attention has been paid to the dependence of the characteristics of individual CID on focal mechanisms. In fact, the majority of large earthquakes occur at deep sea trenches with reverse mechanism. Therefore past studies mostly deal with the ionospheric responses to this particular kind of earthquakes. Here we report observations of the near-field TEC response to three large earthquakes that re-cently occurred in the Kuril Arc with different focal mechanisms, and discuss how their waveforms depend on the coseismic vertical crustal movement patterns.

GPS Data Analysis
Ground-based GPS observations offer a powerful and easy way for remote sensing of the ionosphere. Owing to the dispersive nature of the ionosphere, dual frequency (1.2 and 1.5 GHz) GPS observations provide integral information on the ionosphere by computing the differences in carrier phases recorded by GPS receivers. Methods of TEC calculation have been described in detail in a number of articles (e.g., Calais and Minster, 1995;Afraimovich et al., 2001, and references therein). For convenience, TEC is usually measured in TEC units, TECU (1 TECU=10 16 m −2 ).
Since TEC is an integral parameter, it is impossible to determine the height of TEC disturbance. However, the main contribution to TEC variations would occur around the height of the maximum ionization. This allows us to consider the ionosphere a thin layer located at the height h max of the ionosphere F 2 layer, and TEC represents a point of intersection of a line-of-sight with the thin layer. We trace propagation of CID by subionospheric point (SIP), a projection of an ionospheric piercing point to the earth's surface. In this paper we assumed h max as 300 km.
Low elevation angles tend to enlarge the horizontal extent of the ionospheric region represented by one measurement. Therefore we here used only data with elevations higher than 20 degrees. To eliminate regular diurnal variations of the ionosphere, as well as trends introduced by orbital motion of the satellite, we obtained TEC variations by 1) estimating best-fit polynomials with degree up to seven, and 2) subtracting these polynomials. This works as a high-pass filter to extract variations with periods less than ∼15 min (because we use data period of ∼0.8 hours), which is suitable for studying CID of acoustic wave origin.
The GPS data used in this study all come from the Japanese GPS Earth Observation Network (GEONET, http://terras.gsi.go.jp), operated by Geographical Survey Institute, Japan. It started operation in 1994. Number of GPS stations was ∼100 at the time of the 1994 earthquake, but exceeded 1000 before the 2006/2007 earthquakes. Sampling interval is 30 seconds for all data analyzed here.

General Information on the Earthquakes
The Kuril Islands and Japan are known as one of the regions of the highest seismic activity in the world. Around 20 percent of the world's earthquakes occur in this region. Most large earthquakes in the Kuril and the Northeast Japan arcs occur as interplate low-angle thrust events with magnitudes 7 to 8 that rupture one or multiple asperities (Yamanaka and Kikuchi, 2004). In addition to them, large earthquakes of other mechanisms occur occasionally, e.g. within the subducting slab and at the outer rise of the trench. Here we briefly summarize general characteristics of the three earthquakes studied here based on information available on the web at www.earthquake.usgs.gov and www.eri.u-tokyo.ac.jp/sanchu/SeismoNote/. Coseismic vertical movement patterns, shown in Fig. 5, are all calculated using fault parameters available there assuming an elastic half space (Okada, 1992).
The first earthquake occurred at 13:23 UT (M w 8.3) on October 4, 1994, near Shikotan Island, in the subducting Pacific Plate slab as a high-angle (∼50 deg) reverse faulting resulting in both surface uplift and subsidence of comparable size (Figs. 1, 5(a)). At Shikotan Island, there was approximately 60 cm of subsidence (Kikuchi and Kanamori, 1995). There are two possibilities for the fault geometries that cannot be discriminated with available seismological and geodetic data (Tsuji et al., 1995). However, this ambiguity does not influence the discussion of the present study because the both fault models predict similar vertical coseismic movement patterns characterized by comparable amount of uplift and subsidence.
The second earthquake, an event at 11:14 UT, November 15, 2006, was the largest interplate earthquake (M w 8.2) in the central Kuril Islands since the early 20th century. It occurred in a low-angle thrust fault at the Kuril Trench (Figs. 2, 5(b)), where the Pacific Plate subducts beneath the North American (or Okhotsk) Plate with a velocity of ∼90 mm/year (DeMets et al., 1990;Heki et al., 1999). The main shock occurred at shallow depth within ∼80 km from the trench axis. The vertical crustal movement includes small amount of subsidence, but is dominated by uplift of a meter or more ( Fig. 5(b)).
The third earthquake occurred at 04:23 UT on January 13, 2007 (M w 8.2), in a normal fault within the Pacific Plate (Figs. 3-4, 5(c)), approximately 95 km east-southeast of the November 2006 earthquake. This is a typical outer rise earthquake and is the third largest of this type over the last 100 years, and constitutes a doublet together with the 2006 earthquake (Ammon et al., 2008). The corresponding coseismic crustal deformation is characterized by subsidence of a few meters as shown in Fig. 5(c). Hereafter, we refer to these three earthquakes just as the 1994, 2006, and 2007 earthquakes.

TEC Response to the Earthquakes
First we examine TEC responses to the 2006 earthquake (low-angle reverse faulting at the plate interface), because this is the "typical" large event in subduction zones. About 15 min after the earthquake TEC response appears in records of satellite 20 as TEC fluctuations with amplitude of 0.4-1.0 TECU and period about 500-600 s (Fig. 2). The amplitude of waves differs depending on the azimuth of SIP (compare the right and the left panels), but the response represents a typical compression-rarefaction wave (N -type fluctuations) as a response to propagating shock-acoustic waves (SAW). This agrees with examples in past literatures (e.g., Minster, 1995, 1998;Afraimovich et al., 2001Afraimovich et al., , 2006Heki and Ping, 2005;Astafyeva and Afraimovich, 2006). TEC variations after the 1994 earthquake (high angle reverse faulting within the subducting slab) are presented in Fig. 1. The geometry of satellite 20 allowed us to observe TEC response from two sides (from northwest and southwest) of the epicenter. The shape of the response depends on the azimuth of SIP relative to the epicenter: typical N -type response on the southwestern side (sites 0016, 0024, 0026, 0027, 0031) and "inverted" N -type wave on  the northwestern side (sites 0004, 0005, 0006, 0010). Comparison with Fig. 5(a) shows that we observe the inverted N -type wave at the subsiding side, suggesting that such difference stems from the polarity (i.e. uplift or subsidence) of vertical coseismic crustal movements. Note that the first positive phases of N -like signals at 0016, 0024, 0026, which are in the "neutral zones" between the two sides, are less pronounced than those on the southwestern side of the epicenter (0027, 0031). The TEC responses to the 2007 earthquake (normal faulting at the outer rise), where the coseismic vertical movements were dominated by subsidence (Fig. 5(c)), were rather complicated. The coseismic TEC changes recorded by satellite 1 showed a compound signal as a mixture of different components (Fig. 3(b)): a predominant positive perturbation is followed by an inverted N -wave that diminished after ∼200-300 km. The apparent velocity of the positive disturbance is about 3 km/s (see travel-time diagram in Fig. 3(c)), ∼3 times as fast as the known velocity of propagation of SAW ∼1 km/s (Afraimovich et al., 2001Heki and Ping, 2005). Most likely, this disturbance was caused by the Rayleigh waves. On the other hand, the apparent velocity of the negative disturbance is ∼1.1 km/s (Fig. 3(c)). Figure 4 clearly shows that the disturbances observed with satellite 11 start with negative changes (there the Rayleigh wave signatures are somehow not seen). Similar conclusions can be made, that is, TEC response has a shape of inverted N -wave as we observed above the area of subsidence after the 1994 earthquake, and propagates as fast as ∼1.0 km/s. Thus, our observations suggest that ground subsidence induces ionosphere disturbances starting with negative changes.

Discussions
Using GEONET data, we studied near-field ionospheric responses to the three large (M w > 8) earthquakes in the Kuril Islands, and showed that their waveforms depend on the focal mechanisms. The 2006 thrust earthquake showed a typical compression-rarefaction (N -type) wave, a response to SAW propagating from the uplifted surface of the Earth. TEC response to the 2007 normalfault earthquake was an inverted N -type wave, although it was mixed with another faster-propagating component of Rayleigh wave origin. TEC response to the 1994 high angle reverse-fault earthquake showed both types of N -waves, apparently originating from the uplifted and subsided parts of the Earth's surface.
Generally speaking, the form of the ionospheric response might not exactly follow the form of the source acoustic waves coming from below, since the phases of the source waves along with the way they superpose play important roles in the process of formation of the wave front. The final form of a wave and its amplitude depend also on the geometry of line-of-sight, disturbance wavefront, and especially on the direction of the geomagnetic field (Afraimovich et al., 2001;Heki and Ping, 2005).
It is known that plasma cannot propagate perpendicular to the geomagnetic field lines because of the Lorentz Force exerted on charged particles. This means that in the northern hemisphere the magnetic field prohibits northward propagation of CID (Heki and Ping, 2005). Directions of the geomagnetic field near the epicenters of the three earthquakes are ∼57-58 • in inclination and ∼6-8 • in declination (http://www.ngdc.noaa.gov/geomag/magfield.shtml).
In the case of 1994 earthquake we did observe both northward and southward propagation of CID. However, the northern SIP were located relatively close to the epicenter (just ∼50-150 km away). We hence consider that the wave vector of the CID had not become perpendicular to the geomagnetic field yet when GPS receivers caught these signals. Figure 1 shows that the TEC variations recorded at 0001 and 0004 located farther on the north from the epicenter were much fainter than nearer stations (0005, 0006, and 0010). On the other hand, in the case of southward propagation, we can observe CID farther away from the epicenter with larger amplitudes. We also can see that the westward propagating CID decay faster than those propagating southward in all the earthquakes (Figs. 1-4). These points show that the present cases are consistent with the directivity as found in the 2003 Tokachi earthquake (Heki and Ping, 2005). Waveforms of acoustic perturbations generated at the ground also change during their travel to the ionosphere due to factors including the non-linearity, heterogeneity, and dissipation. If acoustic waves have finite amplitude, the compression part would propagate faster than the rarefaction part because of its higher temperature (i.e. faster sound velocity). This causes steepening of the initial part of the wave, and the whole waveform approaches that of a shock wave (Naugolnykh and Ostrovsky, 1998, section 3.5). Magnetic field may weaken these processes above ∼120 km (Ostrovsky, 2008), and the slope of the wave front decreases and the wavelength extends up to tens of hundreds of kilometers at the height of the F-layer. Atmospheric viscosity also changes the waveform by isolating the wave with frequencies within a window of 3-10 mHz (Blanc, 1985).
In the present study we reported observations of acoustic waves lead by rarefaction parts possibly generated by coseismic ground subsidence. Such waves may not be stable enough to propagate a long distance; because the fasterpropagating compression part may eventually catch up with the slower rarefaction part, they might eventually cancel each other. In fact, our observations show that the amplitudes of the inverted N -type wave are smaller than the normal ones, possibly reflecting their inherent instability. Nonetheless, the present case showed that such inverted Ntype variations still can reach ionosphere in spite of their instability. Although there remain several physical problems in the formation and propagation of the inverted Ntype waves, the observations suggest that the polarity of ground movement (i.e. uplift or subsidence) may cause difference in the waveform. Heki et al. (2006), in studying the near-field ionospheric disturbances by the 2004 Sumatra earthquake, demonstrated that CID has a potential of providing information on focal processes, such as rupture speed and relative magnitudes of multiple asperities. The present study suggests its further potential, i.e. we could use CID to infer focal mechanisms like we have been doing with seismic waves in the solid Earth.