The Persistent Ionospheric Responses Over Japan After the Impact of the 2011 Tohoku Earthquake

In this study, we report the persistent impacts of the 2011 Tohoku earthquake/tsunami on the ionosphere using the ground‐based Global Navigation Satellite System and FORMOSAT‐3/COSMIC total electron content. Multiple unusual ionospheric phenomena, such as ionospheric irregularities, nighttime medium‐scale traveling ionospheric disturbances (MSTIDs), and planar traveling ionospheric disturbances (TIDs), were observed after the emergence of tsunami‐induced concentric gravity waves. The ionospheric irregularities initially developed over the Hokkaido region following the interference of gravity waves at ~8:45 UT. Remarkably, the Perkins‐type nighttime MSTIDs accompanying the planar TIDs were discernible over Japan following the irregularities. By comparing with the tsunami model simulation and ocean buoy observations, it is determined that these planar TIDs, lasting for about 10 hr, were likely related to tsunami ocean waves reflected by seamounts, ridges, islands, and seafloor topography in the Pacific Ocean. Due to the absence of sporadic E layers, we suggest that the coupling between the tsunami‐generated gravity waves and the Perkins instability plays an essential role in initiating the equinoctial nighttime MSTIDs. The long‐lasting tsunami can continuously impact the ionosphere, affecting the nighttime ionospheric electrodynamics and making the conditions conducive for the development of midlatitude nighttime ionospheric irregularities and instabilities.


Introduction
Under quiet geomagnetic conditions, atmospheric wave forcing from the lower atmosphere has considerable impacts on the middle and upper atmosphere (Rishbeth & Mendillo, 2001). Atmospheric waves, such as gravity waves, tides, and planetary waves, generated from different lower atmospheric processes, exhibit various spatial and temporal scales that influence the ionosphere through wave perturbation or electrodynamical processes (e.g., C. H. Lin et al., 2007;H.-L. Liu, 2016). Gravity waves are ubiquitous. They can influence the global circulation as well as background winds, temperature distributions, and constituent structure of the atmosphere, further impacting the global climate and weather forecast. They are perhaps the most crucial dynamical features that can be generated by natural or anthropogenic events on Earth (e.g., Azeem et al., 2015;, Chou, Lin, Shen, et al., 2018C. C. H. Lin, Shen, et al., 2017;. Ionospheric responses to the natural disasters of earthquake/tsunami have been studied for many decades (Bolt, 1964;Leonard & Barnes, 1965;Yuen et al., 1969, and references therein). Transient seismic and tsunami ocean waves can excite acoustic and gravity waves that propagate upward from the Earth's surface to the upper atmosphere, altering the ionospheric plasma, and causing seismotraveling ionospheric disturbances (STIDs). Bolt (1964) first observed two distinct air pressure wave packets associated with seismic and tsunami waves on a Berkeley barogram following the great Alaska earthquake in 1964. The 1964 Alaska earthquake also induced remarkable STID features with periods of approximately 30 s and 90 min in the ionosonde and high-frequency Doppler sounding observations, suggesting that transient perturbations can occur in the ionosphere after large earthquakes (Leonard & Barnes, 1965;Davies & Baker, 1965). Row (1966) further indicated that these short-and long-period oscillations are likely the manifestation of acoustic and gravity waves in the ionosphere. After these initial studies, the scientific community has paid great attention to the vertical coupling between earthquakes and the ionosphere. New efforts have sprung up to provide more observational evidence of earthquake/tsunami-ionosphere coupling (e.g., Afraimovich et al., 2001;Artru et al., 2004;Calais & Minster, 1995;Ducic et al., 2003;Heki & Ping, 2005; J. Y. Otsuka et al., 2006;Rolland et al., 2010;Tanaka et al., 1984). Particular attention has recently been paid to the destructive earthquake/tsunami that occurred off the east coast of Honshu, Tohoku area, Japan, on 11 March 2011. The progress of ground-and space-based instrumentation facilitates further understanding of how the Tohoku earthquake and tsunami affected the upper atmosphere and space environment. The M9.0 Tohoku earthquake excited a devastating tsunami as well as unprecedented ionospheric perturbations. Ground-based Global Navigation Satellite System (GNSS) total electron content (TEC) observations, as well as a cluster of distributed ionosondes in Japan and Taiwan simultaneously recorded remarkable STIDs and concentric TIDs (CTIDs) related to Rayleigh waves and tsunami in unprecedented detail (e.g., Astafyeva et al., 2011;Chen et al., 2011;Maruyama et al., 2011;Rolland et al., 2011;Tsai et al., 2011;Tsugawa et al., 2011). High-frequency Doppler sounders in Taiwan and the Czech Republic also recorded the long-propagating STIDs over thousands of kilometers away from the epicenter (Chum et al., 2012;. Chum et al. (2012) suggested the acoustic wave packets in Doppler shifts are related to seismic waves from the Tohoku earthquake. The CTIDs traveling in a radial direction from the epicenter are the manifestation of tsunami-excited concentric gravity waves in the ionosphere (e.g., Galvan et al., 2012). Multiple satellites and ground-based instruments also captured the tsunami-induced gravity waves propagating throughout the Pacific Ocean and beyond (Azeem et al., 2017;Coïsson et al., 2015;Crowley et al., 2016;Garcia et al., 2013;J. Y. Liu, Chen, et al., 2019;Makela et al., 2011;Smith et al., 2015;Yang et al., 2014). Kakinami et al. (2012) found that the downwelling of the sea surface at the tsunami source area would yield an ionospheric plasma hole due to the plasma downwelling to a lower altitude where the recombination rate is higher. Astafyeva et al. (2013), however, argued that the ionospheric plasma hole was caused by earthquake-generated shock acoustic waves.
Numerical simulations have also been conducted to address the complex processes of earthquake/tsunami-induced acoustic and gravity waves forcing in the middle and upper atmosphere (e.g., Matsumura et al., 2011;Meng et al., 2015;Zettergren et al., 2017). Matsumura et al. (2012) suggested the ground motion could excite acoustic waves trapped in a thermal duct between 50-and 100-km altitudes, forming acoustic resonance extending to the ionosphere with frequencies of 3.7-3.8 mHz as seen in GNSS TEC observations (A. Saito et al., 2011). Zettergren et al. (2017) demonstrated that the nonlinear response by acoustic waves could ultimately lead to the formation of ionospheric TEC depletions as well, which is consistent with the GNSS TEC observations reported by Kakinami et al. (2012).
The studies mentioned above mainly focus on understanding the initial response of the ionosphere to the Tohoku earthquake/tsunami and investigating the characteristics, generation, dissipation, and propagation of STIDs and CTIDs. In this paper, we present the first observation of ionospheric perturbations related to reflected tsunami ocean waves. We outline the effects of the reflected tsunami on the ionospheric electrodynamics and elucidate how tsunami-induced gravity waves leave imprints on the midlatitude nighttime ionosphere. The nonlinear interaction of gravity waves has destabilizing effects on the nighttime ionosphere. Longtime modulation of the tsunami-induced gravity waves in the ionosphere emerges as an essential contributor to modify the ionospheric electrodynamics, leading to the formation of midlatitude nighttime ionospheric irregularities and medium-scale TIDs (MSTIDs).

Ground-Based GNSS TEC
In this study, we use 30-s sampling GNSS observational data in Receiver Independent Exchange format from the GPS Earth Observation Network in Japan. The GPS Earth Observation Network is a compact system of ground-based GNSS receivers that primarily operated for precise point positioning for navigation, surveying, and related activities. It is the largest GNSS network of the world, consisting of 1,200 GNSS permanent stations, providing an opportunity to study ionospheric TEC perturbations with very high spatial and temporal resolution.
We derive the slant TEC from the pseudoranges and phases of the dual-frequency GNSS signals based on the method reported by J. Y. Liu et al. (1996) and determine the instrumental biases by applying an ionospheric spherical harmonic function and least squares method, assuming that the hourly TEC variation is uniform by a given GNSS receiver (Schaer, 1999). The slant TEC is converted to vertical TEC by using a modified single layer model mapping function. The modified single layer model assumes that the ionosphere is a thin, uniform-density spherical shell at the height of 300 km, near the mean altitude of peak electron density. The cutoff elevation angle is set to 20°to avoid multipath errors. Then a fifth-order Butterworth band-pass filter with cutoff periods of 4-30 min is used to extract the characteristics of TIDs for individual vertical TEC observations. As discussed in more detail in section 3, the cutoff periods of 4-30 min were determined based on analysis of the Hilbert spectra.
Butterworth band-pass filter has been widely used to study the ionospheric perturbations related to the Tohoku earthquake (e.g., Chen et al., 2011;Komjathy et al., 2012;. Using this technique allows us to extract wave perturbations within a specific range of wave periods and effectively eliminate the variations due to daily ionospheric background conditions. More distinct wave patterns especially for the small-scale, high-frequency, and weak-amplitude TEC fluctuations can be obtained (e.g., Chou, Lin, Shen, et al., 2018), compared to other detrending methods, such as subtracting the running mean values (e.g., Mrak et al., 2018;Tsugawa et al., 2011).

FORMOSAT-3/COSMIC
The F3/C calibrated radio occultation (RO) TECs (calTECs), another type of observations contained in the level2 ionPrf files in the COSMIC Data Analysis and Archive Center (http://cdaac-www.cosmic.ucar.edu) database, are used to study the ionospheric perturbations associated with the Tohoku tsunami. The F3/C-derived TECs are calibrated to obtain the portion of TEC below the low Earth orbit satellite by subtracting the TECs with positive elevation angles from TECs with negative elevation angles (i.e., the portion of TEC above the low Earth orbit satellite is removed). The calTECs are usually converted to ionospheric electron density profiles as a function of the radius from the Earth's center using Abel inversion (Schreiner et al., 1999). We use the calTECs to study the ionospheric perturbations instead of electron density profiles because the electron density profiles have significant retrieval errors in the bottom side ionosphere due to the spherical symmetry assumption of Abel inversion. Such an assumption will result in artificial plasma caves underneath the equatorial ionization anomaly region (e.g., Pedatella et al., 2015).
Five F3/C satellites were operational during March 2011, giving us about 1,200 successful ionospheric RO profiles per day. F3/C observations provide an opportunity to resolve the vertical perturbations of atmospheric waves in the ionosphere (e.g., Watson & Pedatella, 2018). To extract the vertical disturbances of TIDs, we first interpolated the calTECs into a 1-km vertical resolution grid in tangent point altitude. Then a fifth-order Butterworth filter with a band pass of 20-100 km is applied to remove the background trends.

Hilbert Spectra Analysis
proven method for analyzing nonlinear and non-stationary data and especially useful for estimating low-frequency oscillations. Figure 1 presents the Hilbert spectra of the GNSS TEC on 10-12 March 2011. The Hilbert spectra show significant differences among the three days. The spectra on 11 March show more pronounced high-frequency signals indicated by red rectangles with frequencies ranging from 0.5-5 mHz at 05:50-08:00 UT (14:50-17:00 JST) and~0.6-3 mHz at~12:00-19:00 UT (21:00-04:00 JST). These are notably absent on 10 and 12 March. The initial enhancements from~5:50-08:00 UT are associated with the Tohoku earthquake/tsunami-induced acoustic and gravity waves. In this study, we focus on the second enhancements, which occur between~12:00 and 19:00 UT, and, as will be discussed later, are related to the plasma irregularities, nighttime MSTIDs, and planar TIDs. Based on the HHT spectra, a fifth-order Butterworth band-pass filter with cutoff periods of 4-30 min is applied to extract the wave perturbations in TEC.
There was a geomagnetic storm driven by a corotating interaction region (CIR) on 11 March. The CIR storms are weak in intensity and mainly caused by the high-speed solar winds emanating from a coronal hole. Geomagnetic disturbances at high latitude can be an essential source of large-scale gravity waves and TIDs (LSTIDs) traveling equatorward from the auroral zone (e.g., Richmond, 1978). However, we can rule out the possibility of CIR storm-generated medium-scale and high-frequency waves on 11 March because the small-or medium-scale waves with a shorter period and slower velocity tend to dissipate at high latitude by viscosity, thermal diffusion, and Joule dissipation (cf. Richmond, 1978). Zakharenkova et al. (2016) and  Cherniak and Zakharenkova (2018) revealed that the geomagnetic storm-generated MSTIDs occurred at high latitude and moved slowly as far as of 60-65°N. Most of storm-driven LSTIDs have horizontal phase velocities exceeding 300 m/s and longer wavelength (e.g., Shiokawa, Otsuka, Ogawa, et al., 2002. This demonstrates that only large-scale aurorally generated TIDs can propagate long distances and reach the equatorial region. As the GNSS TEC observations show short-period oscillations, it can be ruled out that these are due to the CIR storm-generated gravity waves at high latitudes.

Observations and Discussions
We present three unusual space weather phenomena, including the ionospheric irregularities, nighttime MSTIDs, and westward planar TIDs several hours after the main Tohoku earthquake/tsunami over Japan, which may have triggered the MSTIDs. These phenomena appear in simultaneous observations of ground-based GNSS TEC and F3/C satellites. The features presented here are unique in that we show the impacts of the reflected tsunami ocean waves on the midlatitude nighttime ionosphere. To better interpret these observations, we divide these space weather phenomena into three sections (sections 4.1-4.3) and address these phenomena separately. Figure 2 shows the time sequence of filtered TEC maps, which illustrates the spatial and temporal evolution of ionospheric disturbances after the Tohoku earthquake. Figure 2a shows the Tohoku earthquake-excited distinct CTIDs emanating outward over Japan at 07:45 UT (16:45 JST) (e.g., J. Y. Tsugawa et al., 2011). These CTIDs, with horizontal phase velocities of 200-310 m/s and wavelengths of 200-300 km, are the manifestation of tsunami-excited concentric gravity waves in the ionosphere (e.g., Galvan et al., 2012;Tsugawa et al., 2011). Interestingly, we observed clear planar TIDs emanating eastward/southeastward from the Pacific coast of Japan during~07:15-08:30 UT (16:15-17:30 JST). The eastward planar TIDs aligned with the Pacific coast of Japan could be related to the tsunami backscattered by the Japanese mainland and excited eastward propagating planar gravity waves. Note that the eastward planar TIDs have a horizontal phase velocity of~200 m/s, horizontal wavelength of~240 km, and period of 20 min, which are identical to the characteristics of CTIDs.
The propagation of gravity waves can cause an inhomogeneity in ionospheric conductivity. The inhomogeneous conductivity within the background wind field will be electrically polarized to satisfy the quasineutrality condition (i.e., ∇ · J = 0), which can further induce plasma irregularity through E × B drifts. Additionally, gravity wave wind perturbations δU are able to drive these small-scale ionospheric irregularities through dynamo currents δJ = Σ F p (δU × B) (e.g., Hysell et al., 2018;Tsunoda, 2010). The dynamo currents flow across the geomagnetic field lines where the Pedersen conductivity is inhomogeneous. The resulting polarization electric fields (δE) can trigger ionospheric irregularities.
The gravity wave-driven δE in the F layer should depend on the ratio of field line-integrated Pedersen conductivity in the E and F regions as given by (Abdu et al., 2015) δE¼ where Σ F p and Σ E p are field line-integrated Pedersen conductivities in the E and F regions. The conductivity ratio will increase toward sunset due to the decrease of E region conductivity, leading to an increase of δE, which can drive plasma irregularities. The plasma irregularities cannot develop during the daytime because the ionospheric E region Pedersen conductivity is large, meaning that the induced δE will short out. However, the induced δE is sustainable without shorting out at sunset. This process resembles the gradient drift (E × B) instability (Kelley, 2009) and the Perkins instability (Perkins, 1973) in the midlatitude ionosphere, which could be a key factor for midlatitude spread F (e.g., C. A. Miller et al., 1997). Coster et al. (2017) observed wave-like features of enhanced TEC propagating eastward over the Rocky Mountains during the solar eclipse over North America on 21 August 2017, suggesting that large temperature and pressure changes could locally produce atmospheric waves. The plasma irregularities also developed simultaneously near the same region (Pradipta et al., 2018). Pradipta et al. (2018) suggested that the development of plasma irregularities may be related to Rayleigh-Taylor instability or E × B instability or a momentary surge in field-aligned current due to the reduction of E and F region conductivities. The coincident observations demonstrate that the interaction of atmospheric waves and the reduction of the E region Pedersen conductivity creates a favorable condition for the development of plasma irregularities.
Note that the plasma irregularities extended to a broader patch, and the southeastward TIDs ( Figure 2d) emanating from the irregularity patch was discernible around 10:17 UT (19:17 JST), implying that the localized plasma irregularities could be a potential source to excite TIDs. A similar phenomenon was also observed by Chou, Lin, Yue, Chang, et al. (2017), who identified small-scale TEC striations following the typhoon-excited concentric gravity waves in the ionosphere. The MSTIDs subsequently developed from the small-scale TEC striations. Pradipta et al. (2018) also observed TID signatures following the midlatitude spread F from the same ionosonde station during the total eclipse over North America on 21 August 2017. The underlying physics for plasma irregularities to generate TIDs remains a puzzle. There are two plausible scenarios to explain the formation of TIDs. The first scenario is that the plasma irregularities could excite the TIDs due to localized heating. This may occur due to Joule heating effects by the small-scale electric field fluctuations within the irregularities (e.g., Codrescu et al., 2000;Kagan & Kelley, 2000;Kelley, 2009). Polarization electric fields within the midlatitude ionospheric instabilities could be around 1-5 mV/m (e.g., Hysell et al., 2018;. Kagan and Kelley (2000) and Kelley (2009)  proposed a thermal instability, suggesting that electric field perturbations caused by either a neutral wind or external electric field could lead to plasma heating and instability while the background dynamo current and polarization electric field are in the same direction (see Figure 6.44 of Kelley, 2009). The Joule heating can reduce the plasma recombination rate (cf. Schlegel, 1982), which in turn increases ionospheric plasma and ion-neutral collision frequency, producing more Joule heating. The local heating can cause temperature and neutral winds perturbations, leading to ionospheric disturbances (cf. C. C. H. Lin, Shen, et al., 2017. An alternative possibility is that gravity waves generated by the tsunami may be a potential source to drive the southeastward TIDs. Gravity waves breaking or dissipating can also cause heating/cooling in the thermosphere, creating body forces to excite secondary gravity waves in the ionosphere (e.g., Vadas, 2013;Vadas & Crowley, 2010).

Equinoctial Nighttime MSTIDs
Following the plasma irregularities, Figure 3 shows that distinctive band structures with northwest-southeast wavefront alignment, referred to as nighttime MSTIDs, also occurred off the east coast of Hokkaido (Figure 3a) at~11:50 UT (20:50 JST). The MSTIDs initially have a peak amplitude of 0.34 TECU (~2% of background TECU) and propagate southwestward with increasing amplitudes of 0.6-0.9 TECU (~4-6% of background TECU). Meanwhile, the plasma irregularities gradually become evanescent. These wavy structures eventually reached~25°N and then faded away at~16:30 UT (01:30 JST) (not shown). Notably, the planar TIDs were visible over Japan during the emergence of MSTIDs. The planar TIDs with a relatively weak amplitude of~0.15 TECU (~1% of background TECU) propagating westward, northwestward, and southwestward are likely related to reflected tsunami-triggered gravity waves. The details of the planar TIDs will be discussed in the next section.

Space Weather
Nighttime MSTIDs usually have periods ranging from several minutes to less than an hour, velocities of 50-300 m/s, and amplitudes of up to few TECU (e.g., Hernández-Pajares et al., 2006). We estimated the morphological characteristics of MSTIDs using the keogram as a function of the time and distances by organizing the filtered TECs along the solid line shown in Figure 3d. Figure 4 shows the keogram that the MSTIDs have a horizontal phase velocity of~170 m/s, a period of~26 min, and a horizontal wavelength of~265 km, which agrees well with previous statistical studies of MSTIDs over Japan (e.g., Tsugawa, Kotake, et al., 2006). Note that the observed MSTIDs have higher amplitudes than the planar TIDs. Tsugawa, Kotake, et al. (2006) reported that the nighttime MSTIDs have higher amplitudes (1-5%) than the daytime MSTIDs (1-2%), suggesting the electrodynamic force of the Perkins instability is critical. The fact that the equinoctial nighttime MSTIDs fade away around~25°N is also consistent with previous studies (e.g., Shiokawa, Otsuka, Ejiri, et al., 2002). Shiokawa, Otsuka, Ejiri, et al. (2002) and Narayanan et al. (2014) suggested the latitudinal constraint of MSTIDs is associated with the enhancement of field line integrated Pedersen conductivity in the equatorial ionization anomaly crests that reduce the polarization electric fields within the MSTIDs.
The Perkins instability was believed to be the most likely mechanism to explain the specific wavefront alignment of MSTIDs (Perkins, 1973). However, a coupled electrodynamical process is usually necessary to reinforce the instability due to its small linear growth rate (e.g., Cosgrove & Tsunoda, 2004;Kelley & Fukao, 1991).  observed positive and negative charge accumulation at the edge of the most intense airglow bands. It is suggested that a southeastward wind could drive a northeastward Pedersen current flowing along regions where inhomogeneity of Pedersen conductivity, developing polarization electric fields at the edge of the band structures, triggering the Perkins instability (e.g., Otsuka et al., 2007). The polarization electric field could be perpendicular to the wavefront of the MSTIDs. The Horizontal Wind Model-14 (Drob et al., 2015) exhibits southeastward wind patterns over Japan during these periods, suggesting that the Perkins instability may account for the observations. The maximum linear growth rate (γ) of the Perkins instability can be expressed as (Perkins, 1973) γ where E 0 , B, H n , I, θ, g, and ν in are the background electric field, geomagnetic field, neutral scale height, magnetic inclination angle, the angle of background electric field from the geomagnetic east, and plasma density-weighted height-integrated ion-neutral collision frequency. Nighttime MSTIDs have a high occurrence rate during solstices and low solar flux years (e.g., Kotake et al., 2006;. A potential explanation is that the neutral density in the upper atmosphere is small during solstices and low solar flux years, which could result in lower ion-neutral collision frequency. The linear growth rate is inversely proportional to the plasma density-weighted height-integrated ion-neutral collision frequency. This fact supports the double peaks of MSTID occurrence in June and December solstices. The linear growth rate in equinox is usually smaller. Therefore, an additional driver is required to accelerate the growth of the equinox MSTID. A CIR storm occurring during these periods can increase the linear growth rate by generating penetration electric field to amplify the background electric field. However, Nishioka et al. (2009) reported a small linear growth rate (1.6 × 10 −4 s −1 ) of super MSTIDs during a super geomagnetic storm (Dst < −250 nT), demonstrating that the CIR storm itself might be insufficient to drive the MSTIDs with the amplitude higher than usual.
A coupled electrodynamical system between the sporadic E (Es) layer and F region is the predominant mechanism to explain the faster growth of nighttime MSTIDs (e.g., Cosgrove & Tsunoda, 2004;Tsunoda, 2006;Yokoyama et al., 2009). Otsuka et al. (2008) investigated the relationship between MSTID and Es layer activities in summer nights. The Es activity is closely correlated with the MSTID activity, suggesting that the occurrence of Es is critical for the generation of MSTIDs. The higher magnitude of the Es layer in solstices observed by F3/C also supports the increased occurrence of MSTID in solstices (Yeh et al., 2012). In this case, a persistent high conductivity E layer is required to grow the polarization electric fields between the E and F regions (Otsuka et al., 2007). We examine the Es layer by using ionograms from the National Institute of Information and Communications Technology in Japan. Since the Es layer may extend over a wide range (cf. Yeh et al., 2012Yeh et al., , 2014, the F3/C excess phase (L1 and L2 channel) and signal-to-noise ratio of the L1 (L1 SNR) from atmPhs files are also used to examine these features (cf. Yue et al., 2015). Figure 5 shows the ionograms at Kokubuji and Yamagawa stations (top) and the F3/C RO normalized L1 SNR and excess phase (bottom). The Es features are nearly absent from both observations. Only the fourth F3/C RO event over the Chinese region shows more significant fluctuations in L1 SNR and excess phase in comparison with background fluctuations, suggesting the occurrence of an Es layer. Both observations seemingly rule out the possibility that the Es-F layer coupling is responsible for the observed MSTIDs.

Space Weather
Consequently, due to the absence of Es layer, gravity waves are considered as a potential driver of the nighttime MSTIDs by coupling with the Perkins instability (e.g., Kelley & Fukao, 1991;C. S. Huang et al., 1994;C. A. Miller et al., 1997;Kelley, 2011;Chou, Lin, Yue, Chang, et al., 2017. The southwestward TIDs from the east coast of Hokkaido could be related to the local heating or gravity waves excited by reflected tsunami ocean waves. The southwestward TIDs modulate with the tsunami-induced TIDs, subsequently coupling with the Perkins instability to generate the equinox nighttime MSTIDs with increasing amplitude (Figures 3f-h) (cf. Figure 5 of Chou, Lin, Yue, Chang, et al., 2017). The westward and southwestward propagating TIDs can more efficiently promote the growth rate of the Perkins instability (Chou, Lin, Yue, Chang, et al., 2017;, making the rare observation of equinoctial nighttime MSTIDs. Notably, the MSTIDs have a dominant horizontal wavelength of~265 km that is comparable to the tsunami-induced gravity waves of~290 km, implying the tsunami-induced gravity waves are critical.  performed model simulations of MSTIDs by considering gravity waves and an oscillating electric field. They suggest that the oscillating electric fields are critical in the generation of MSTIDs.  further showed the geomagnetic conjugacy of nighttime MSTIDs to support the electrodynamical effect. Simulations of gravity waves through wind perturbations lead to a strong conjugate effect in the ionosphere as well (Wu et al., 2015), suggesting that, in addition to the Es-F coupling, gravity wave seeding can also explain the generation of nighttime MSTIDs.
It is noteworthy that the mechanism behind the Es-F layer coupling is not well established. Cosgrove and Tsunoda (2004) and Yokoyama et al. (2009) suggested the polarization electric fields of Es layer instability could map to the F layer, considerably enhancing the growth rate of the Perkins instability. On the other hand, S. Saito et al. (2007) argued that the horizontal scale of quasiperiodic echoes (~10 km) in the E region is much smaller than that of MSTIDs (~150 km), which is different from the postulated situation in the previous model simulations of Es-F layer coupling. Future investigation further necessitates knowledge of understanding how the scale size in both E and F regions can reconcile with each other. To understand the origins of the planar TIDs, we further compare them with the Method Of Splitting Tsunami model simulations. Figure 7 shows the tsunami propagation and evolution from Movie S1 of Tang et al. (2012). Multiple tsunami wavefronts are seen to propagate back toward Japan with nearly concentric or planar patterns. These reflected waves are due to the initial tsunami waves encountering seamounts, ridges, islands, and interacting with seafloor topography in the North Pacific Ocean. Figure 7a shows the southwestward (yellow rectangle) and northwestward (red rectangle) propagating tsunami backscattered by the Emperor Seamounts and Mid-Pacific Seamounts at 10:31 UT (19:31 JST). Comparison with Figure 3a shows that these reflected waves have similar wavefront orientations as the northwestward and southwestward planar TIDs. The southwestward tsunami-excited gravity waves could play an important role in seeding the nighttime MSTIDs. Figures 7b and 7c show that the tsunami (red rectangles) backscattered by Mid-Pacific Seamounts at 13:26 UT (22:26 JST) and by Marshall Island at 14:46 UT (23:46 JST), respectively. The tsunami wave patterns are slightly different. The reflected tsunami waves in Figure 7b have a planar wave pattern, which is identical to the planar TIDs shown in Figures 3d and 6a. The reflected waves in Figure 7c display a concentric wave pattern similar to the TIDs wavefronts shown in Figures 6b-6d. Figure 7d shows the tsunami backscattered by South America reaching Japan at 03:56 UT (12:56 JST) on 13 March, two days after the earthquake. The wavefront orientations are also similar to the planar TIDs shown in Figures 6e and 6f. Overall, the planar TIDs reflect the wavefront orientation of the reflected tsunami waves at different times, suggesting that their source is the tsunami. It is noteworthy that the time latency between the tsunami and planar TIDs maps is approximately 1-2 hr, which is reasonable time delay since the tsunami-induced gravity waves may reach the ionosphere~1-2 hr after excitation at the surface (Artru et al., 2005;Vadas et al., 2015).

Reflected Tsunami-Induced TIDs
To further understand the characteristics of the planar TIDs, Figure 8 shows the keogram of the planar TIDs by organizing the filtered TEC along the solid line shown in Figure 6d. The keogram displays long-lasting coherent wave patterns for about 10 hr (10:00-20:00 UT; 19:00-05:00 JST) with a horizontal phase velocity of~230 m/s, horizontal wavelength of~290 km, and period of~21 min. These characteristics are nearly identical to some of the tsunami-generated TIDs (e.g., Azeem et al., 2017;, 2011Makela et al., 2011), however, have some differences in comparison with the CTIDs (Figure 2a). We note that the reflected tsunami-induced TIDs have a nearly coherent wave pattern, while Figure 2a shows the CTIDs are dispersive. The horizontal wavelengths in the inner rings are much smaller than in the outer rings. Such a tendency can be seen in convectively generated gravity waves (e.g., Azeem et al., 2015;. Gravity waves launched from a point source would spread out horizontally and vertically, leading to conically shaped wave phase fronts (J. Yue et al., 2013). Small-scale waves usually have steep propagation angles that can travel shorter horizontal distances (J. Yue et al., 2009). Note that the southwestward propagating CTIDs in the Pacific Ocean are colocated with the tsunami wavefronts (e.g., Galvan et al., 2012) and have similar horizontal wavelengths as the planar TIDs. The small-scale CTIDs may be related to the ground motion of the Tohoku earthquake (e.g., Astafyeva et al., 2011), which could act as point source partially contributing to the dispersive CTIDs. Figure 9a shows eight trajectories of F3/C radio occultation calTEC profiles, recording ionospheric perturbations associated with reflected tsunami ocean waves (#1-#4) and reference observations (#5-#8). The solid and dashed lines shown in Figure 9b indicate the filtered and calibrated TEC corresponding to Figure 9a. The blue lines indicate the observations on 11 March 2011, and the red lines indicate the reference observations with similar occultation geometries. F3/C detected the TEC perturbations with a peak amplitude of 3.3 TECU increasing exponentially from 110 to 300 km and gradually decreasing amplitude above 300 km. The corresponding wavelet spectrograms (Figure 9c) show that the perturbations have vertical wavelengths ranging from~50 to 100 km between 150-and 500-km altitudes. This agrees well with the GOCE and GRACE satellite observations that the tsunami-generated gravity waves can reach 270-to 450-km altitude without severe dissipation (Garcia et al., 2013;Yang et al., 2014). The second profile only shows weak fluctuations, which could be related to the track of tangent points being nearly parallel to the planar wavefronts (e.g., Meng et al., 2015). The fourth profile shows greater amplitudes around~17:29 UT (02:29 JST), which is consistent with the enhanced ground-based TEC perturbations shown in Figure 8. Note the F3/C observations show greater wave amplitude in comparison with the ground-based GNSS

10.1029/2019SW002302
Space Weather TEC observations, which can be attributed to the relative orientation between the ray path geometry of the F3/C occultations and wavefront orientation of gravity waves. Some cancelation effect would occur when the ray path of GNSS intersect different phases within the gravity waves. In particular, the occultations with a ray path of F3/C parallel to the wavefronts of gravity waves usually can obtain higher amplitude.
We can illustrate the planar TIDs by using the Boussinesq dispersion relation for linear gravity wave theory in the isothermal condition; the dispersion relation is given by (Hines, 1960) where m and k are the vertical and horizontal wave number, N is Brunt-Väisälä frequency, ω I = k(c h − u) is intrinsic frequency, and c h and u are the horizontal phase velocity of gravity wave and background wind speed in the direction of wave propagation. We can assume the average Brunt-Väisälä frequency N ¼ 2π 8 min from 0 to 300 km using empirical neutral atmosphere parameters from NRLMSISE-00 (Picone et al., 2002). The Horizontal Wind Model-14 shows that the background winds are southeastward with a velocity of 50 m/s. The wind direction favors the northwestward propagation of reflected tsunami-induced gravity waves. Based on these parameters, the vertical wavelength λ v ¼ 2π m À Á estimated from dispersion relation is 120 km, which is roughly comparable to the F3/C observations. We further make use of observations of water column height (WCH) variations from Deep-ocean Assessment and Reporting of Tsunami (DART) stations of 21413, 21418, and 21419 to investigate the relationship between the tsunami and planar TIDs. Figure 10 shows the locations of the DART stations (Figure 10a), the corresponding WCH variations on two consecutive days (11 and 12 March 2011) (Figures 10b, 10e, and 10h), and the filtered WCH with band-pass filtering of 4-60 min (Figures 10c,  10f, and 10i) on 11 March and the wavelet spectrogram (Figures 10d, 10g, and 10j). The WCH observations on 11 March display significant fluctuations from 10:00-20:00 UT (19:00-05:00 JST) in comparison with the reference day (12 March), indicating the propagation of tsunami ocean waves. We choose these periods because the planar TIDs are discernible after 10:00 UT (19:00 JST) (Figure 8). It should be mentioned that the WCH observations on 11 and 12 March are at 1-and 15-min temporal cadence, respectively. This does not account for the lack of long-period fluctuations on 12 March, and significant long-period tsunami ocean waves are seen in the DART observations sampled at a 15-min temporal cadence on 11 March (not shown).  (Figures 10d and 10g). The time difference of these waves between 21413 and 21418 may be due to the wave propagation. We note that the wave periods are slightly different among these stations, probably because multiple tsunami ocean waves coming from different directions with various periods affect the wave pattern and result in a broad spectrum. These reflected tsunami ocean waves arriving at~13:30 UT could be related to the reflected tsunami waves as shown in Figure 7b. They can be the source of the planar TIDs at~15:00-19:00 UT (24:00-04:00 JST) as shown in Figure 8. Likewise, Figure 10j shows the reflected tsunami waves with similar periods of 15-40 min around 10:00-12:00 UT (19:00-21:00 JST). These reflected tsunami ocean waves propagated southwestward ( Figure 7a) and triggered the southwestward propagating planar TIDs (Figure 3a). There is approximately 1-1.5 hr delay between the TIDs (Figure 8) and WCH perturbations (Figure 10g), which can be due to the propagation of gravity waves from the ocean to the ionosphere. Note that the WCH observations show a long period of background variation on both days (Figures 10b, 10e, and 10h). It could be related to the lunisolar tides, which has an average period of 12h25m (i.e., a half-rotation of the Earth relative to the Moon). Makela et al. (2011) reported that the tsunami-induced ionospheric perturbations precede the tsunami ocean waves by approximately 1 hr. They suggested the early arrival waves in the ionosphere are related to small oscillation tsunami waves generated by prerupture processes. In this study, the small-scale waves were also identified around 11:00-15:00 UT (20:00-24:00 JST). However, it is difficult to disentangle the preceding 10.1029/2019SW002302 Space Weather ionospheric disturbances due to the complicated sea level and ionospheric conditions that have already been disturbed by the main earthquake rupture and the tsunami.
A discrepancy of wave periods between the WCH perturbations (15-40 min) and TIDs (~21 min) exists. It may be due to the tsunami exciting a broad spectrum of gravity waves (Vadas et al., 2015). Short-period waves usually have longer vertical wavelengths that can survive in the ionosphere before dissipation by molecular diffusion (Vadas, 2007). In contrast, long-period tsunami waves have shorter vertical wavelengths that will quickly dissipate before reaching the ionosphere. This also explains the absence of long-period waves in the GNSS TEC observations. Since the gravity waves propagate upward obliquely, we can estimate the horizontal distance R traveled by gravity waves using the following relation (e.g., J. Yue et al., 2009): where z and β indicate the vertical distance from the source to the ionosphere and zenith angle between vertical and wave vector. Assuming z = 300 km from the ocean surface to the ionosphere, N ¼ 2π 8 min , and the R can be estimated to be~726 km. This indicates the gravity waves generated by tsunami could be detected in the ionosphere~726 km away from the source.
We further estimate the vertical group velocity of gravity wave using the following equation: where t represents the propagation time from source to the ionosphere, β ¼ tan − 1 R z ¼ 67:5°, and k ¼ 2π 290 km , then the V g and t can be estimated to be~0.081 km/s and~62 min. The dispersion relation suggests that the tsunami-generated gravity waves may take about 1 hr to propagate from the source to the ionosphere, where they generate TIDs near the east coast of Japan. The time latency roughly agrees with the~1-1.5 hr delay between planar TIDs and WCH perturbations, as well as the GOCE observations of the propagation time of gravity waves from the ocean surface to the GOCE altitude (270 km) (Garcia et al., 2014). The estimated horizontal distance also agrees with the TEC observations that the planar TIDs appeared near the blue dashed circle with a radius of 726 km from 213418 station (Figure 10a), suggesting that the reflected tsunami ocean waves may be responsible for the planar TIDs.  suggested that the small-and medium-scale gravity waves from the lower atmosphere can reach the thermosphere and dissipate at altitudes of 120-250 km. As a result, horizontal forces were created in the thermosphere when gravity waves dissipated momentum and energy. These processes can lead to neutral wind perturbations and generation of secondary gravity waves. The observed TIDs could be a manifestation of secondary gravity waves in the ionosphere (e.g., Nishioka et al., 2013). In this study, the reflected tsunami-triggered TIDs are likely primary gravity waves since the TIDs observations agree with the dispersion relation. Vadas et al. (2015) simulated the ocean wave-excited gravity waves and found that they reach ionospheric heights (~250 km) after~75-80 min. Azeem et al. (2017) suggested that the Tohoku tsunami-induced TIDs over the United States originated from the west of the Pacific coast of the United States by applying a 3-D reverse ray trace model. These studies revealed that the tsunami-induced gravity waves could directly propagate to the ionospheric F region from the ocean surface without significant dissipation.

Summary
In this paper, we investigated the long-lived ionospheric responses after the 2011 Tohoku earthquake and tsunami by using ground-based GNSS and F3/C TEC observations. Multiple ionospheric features such as planar TIDs, plasma irregularities, and nighttime MSTIDs developed following the earthquake-and tsunami-induced CTID. We suggest that the widespread tsunami waves reflected by seafloor topography in the Pacific Ocean could be the primary driver to trigger the planar TIDs through gravity wave perturbations. The modulation of planar TIDs in the ionosphere likely plays a crucial role in seeding the small-scale plasma irregularities and nighttime MSTIDs. The main results of the present work are as follows: 1. Plasma irregularities organized off the east coast of Hokkaido are likely related to the tsunami-induced gravity waves. They are capable of driving these plasma irregularities through dynamo currents, leading to the growth of plasma irregularity through E × B drifts. 2. TIDs originating from the irregularity patch could be related to local plasma heating effects caused by turbulent electric field perturbations within the gravity waves and plasma irregularities. Alternatively, primary or secondary gravity waves related to the tsunami could be another potential source to excite TIDs. The excited TIDs are further modulated by the tsunami-induced southwestward and westward propagating TIDs, and subsequently amplified by the Perkins instability to develop the equinoctial nighttime MSTIDs. 3. Equinoctial nighttime MSTIDs that are rarely seen in equinox propagating southwestward with a horizontal phase velocity of~170 m/s, a period of~26 min, and a horizontal wavelength of~265 km, which are identical to the solstitial nighttime MSTIDs. The peak amplitude varied from~0.34 to~0.9 TECU while propagating southwestward, implying the contribution of the Perkins instability with external driving forces from tsunami effects. 4. Gravity wave seeding is the most plausible mechanism to explain the rare observation of equinoctial nighttime MSTIDs in the absence of an Es layer. The horizontal wavelength of MSTIDs is similar to the tsunami-induced TIDs, suggesting that electrodynamical coupling between the tsunami-generated gravity waves and the Perkins instability is critical for the formation of equinoctial nighttime MSTIDs. 5. The planar TIDs last for more than 10 hr and have a horizontal phase velocity of~230 m/s, horizontal wavelength of~290 km, vertical wavelengths of 50-100 km, and a period of~21 min. They are most likely related to the tsunami ocean waves reflected by seamounts, ridges, islands and seafloor topography in the North Pacific Ocean based on the following reasons: (1) The planar TIDs have similar horizontal wavelength, speed and period as the tsunami-induced TIDs (e.g., Azeem et al., 2017;J. Y. Liu, Tsai, Ma, et al., 2006, 2011Makela et al., 2011); (2) the wavefront orientations of the planar TIDs is similar to the reflected tsunami ocean waves; (3) the~1-1.5 hr delay between the planar TIDs and DART WCH variation agrees with the theoretical propagation time of gravity waves from the ocean to the ionosphere; and (4) the planar TIDs appeared on the east coast of Japan, agreeing well with the theoretical horizontal distance of gravity wave propagation from the DART station to the ionosphere.
This study reports that the ionospheric responses to the Tohoku tsunami and its subsequent disturbances last for many hours, and the induced gravity waves can significantly impact the ionosphere and thermosphere system. The extended effects have important implications for space weather applications and should be taken into consideration for future space weather forecasts. Gravity waves are a potential source to generate ionospheric irregularities such as Es layer instability and equatorial and midlatitude spread F (e.g., Haldoupis, 2012;Hysell et al., 2018;Tsunoda, 2010), which can cause amplitude scintillations and phase fluctuations in transionospheric radio signals, leading to degradation in GNSS performances (e.g., . Significant scintillations typically occur in the equatorial and auroral regions, as the intense ionospheric disturbances are usually limited to low and high latitudes (e.g., Cherniak & Zakharenkova, 2018;Zakharenkova et al., 2019). At midlatitudes, localized ionospheric scintillations occur most frequently during summer nighttime and solar minimum years, in association with the well-known midlatitude spread F (e.g., Hajkowicz & Minakoshi, 2003;Kotake et al., 2006). Yoon and Lee (2014) investigated the potential impact of midlatitude nighttime MSTIDs and irregularities on GPS-based navigation systems. They found loss of lock of GPS signals frequently occurred during the MSTID events, especially when the ionospheric spatial gradient is significant, leading to potential threats to users of ground-and space-based augmentation systems. These studies demonstrate that the observed nighttime MSTIDs and small-scale ionospheric irregularities can be a source of midlatitude ionospheric scintillations. They can also cause large ionospheric spatial gradients, leading to errors in GNSS radio occultation retrievals (cf. Pedatella et al., 2015), which in turn result in errors in space weather forecasts (e.g., Chen et al., 2017;Hsu et al., 2018;C. Y. Lin, Matsuo, et al., 2017;Rajesh et al., 2017).
Additionally, although the observed ionospheric phenomena in this study lead to perturbations of less than 1 TECU, and their effects are not as dramatic as equatorial spread F. Hernández-Pajares et al. (2006) demonstrated that the presence of MSTIDs leads to errors in the ionospheric correction for GNSS receivers. Even small TEC fluctuations (~0.25 TECU) caused by the MSTIDs can still degrade the performance of precise navigation. They further developed an approach for correcting MSTID effects on precise real-time kinematic and tropospheric determination (Hernández-Pajares et al., 2017). Consequently, the impacts of TIDs and plasma irregularities reported on in this study have a nonnegligible effect on GNSS navigation and applications.
Furthermore, the ground-based GNSS networks have proved to be a useful tool for the tsunami early warning system (Blewitt et al., 2006;Song, 2007). Recent studies demonstrated that the origin of a tsunami could be located utilizing ground-based GNSS TEC observations (Kamogawa et al., 2016; J. Y. Liu, Tsai, Ma, et al., 2006;J. Y. Liu, Lin, et al., 2019;Tsai et al., 2011). Rakoto et al. (2018) showed that TEC observations could be inverted to derive tsunami wave heights that are in reasonably good agreement with DART observations. Nowadays, more than thousands of ground-based GNSS stations and space-based GNSS RO satellites have been routinely operating in real time, shedding light on the potential for joint ground-and space-based GNSS observations for tsunami detection and warning, and advancing our understanding of how tsunamis drive space weather effects.