# Preseismic Perturbations and their Inhomogeneity as Computed from Ground- and Space-Based Investigation during the 2016 Fukushima Earthquake

Authors: Sagardweep Biswas; Subrata Kundu; Sudipta Sasmal; Dimitrios Z. Politis; Stelios M. Potirakis; Masashi Hayakawa
Journal: Journal of Sensors
 (2023)

---

## Abstract

We present the atmospheric anomalies instigated through seismogenic sources by a multichannel observation using ground- and satellite-based systems. This study emphasizes the seismic event which happened on the east coast of Japan, near the Fukushima Prefecture on November 21, 2016 (in UTC), with a magnitude of 6.9 and a depth of 11.4 km. We mainly focus on the atmospheric and ionospheric irregularities via acoustic and electromagnetic channels originating from earthquakes in the process of the lithosphere, atmosphere, and ionospheric coupling (LAIC) mechanism. In the acoustic channel, we study the seismogenic atmospheric gravity wave (AGW) which perturbs the local lower atmosphere. The observation of nighttime fluctuations in the very low frequency (VLF) signals and total electron content (TEC) is used to investigate the atmospheric perturbation through the electromagnetic channel. For the ground-based observations, a VLF signal network consisting of 5 receivers in Japan is used to study by recording the VLF amplitude transmitted from the Japanese transmitter JJI (22.2 kHz). VLF nighttime fluctuation is used to check the unusualities due to the earthquake. Preseismic wavelike structures having periods of AGW are observed in the nighttime signal. Direct investigation of such AGWs is done by computing the potential energy related to AGW from the sounding of the atmosphere using broadband emission radiometry (SABER) temperature profiles mounted on the Thermosphere Ionosphere Mesosphere Energetics and Dynamics (TIMED) satellite. Ionospheric TEC inspection is done by using a ground-based global navigation satellite system (GNSS) receiver from the International GNSS Survey (IGS) station MIZU in Japan and observing anomalies in diurnal TEC around 6 and 10 days prior to the earthquake. We also obtain the wavelike structure of AGW from the small-scale fluctuation of TEC using wavelet analysis. All the parameters are found to be preseismic for this earthquake; the acoustics channel gives more consistent outcomes than the electromagnetic channel.


---

## Body

## 1. Introduction

In several investigations in the past few decades, it has been well established that a major disruption in the normal atmospheric processes gets perturbed due to seismic hazards and to investigate this, a hypothesis is proposed, named as lithosphere-atmosphere-ionosphere coupling (LAIC) [1–11]. According to LAIC, the seismogenic impression can be estimated by three channels: (a) the chemical channel, (b) the acoustic channel, and (c) the electromagnetic channel. In the chemical channel, radon plays an active role, as proposed by Pulinets and Ouzounov [11]. In the acoustics channel, intensification in the energy of atmospheric gravity waves (AGWs) is the main observable to understand pre- and coseismic anomalies [12–16]. In the electromagnetic channel, a group of parameters are found to be useful, such as anomalies in the very low frequency (VLF) radio signal [1, 9, 17], unusual emergence of the ultra-low frequency (ULF) signal [18, 19], irregularities in the critical frequency (foF2) [20–22], anomalies in the ionospheric total electron content (TEC) [23–25], and many more. In this manuscript, we use two such channels to investigate the preseismic signatures during the Fukushima earthquake.The electromagnetic impression of preseismic irregularities from the subionospheric VLF radio technique as a ground-based observation became highly convincing during the Kobe earthquake in 1995. Hayakawa et al. [1] reported a sudden shift in the VLF terminator times a few days before the earthquake. After that, extensive statistical studies evidenced the strong correlation between the preseismic mechanism and the VLF signal disturbances [26–28]. The statistical studies had been applied by varying the depth and magnitudes of earthquakes [29–32], which yielded that shallow earthquakes (<40 km), and earthquakes above a magnitude of M ≥5.5 gave a good correlation by this method. A recent study also found the existence of criticality conditions in the lower ionosphere during earthquake activity from VLF observations [33, 34].Conversationally, three basic techniques are found to be useful for the analytical solutions using the VLF signal to detect seismogenic modulation. They are (a) the terminator time method [5, 17], (b) D-layer preparation and disappearance time (DLPT and DLDT) abnormalities [9] and (c) nighttime fluctuation [6, 29]. Among those methods, nighttime fluctuation has the significant advantage of exemption of the superposition of the effects of solar activities in the daytime. The method of this study was based on the nighttime portion of the signal by the standard statistical measurements of abnormal signal behavior from its normal condition [6, 27, 35, 36]. The early method is developed by subtracting the individual nighttime data with the average data to observe the abnormal behavior which refers to a trend. After that, the statistical standard deviations and fluctuations are introduced [28, 31]. Another method proposed by Ray et al. [37, 38] is by estimating the quiet condition with a statistical standard deviation filter and showing a statistical correlation between abnormalities in the signal deviation and earthquake effective magnitudes.Asano et al. [39] found a significant preseismic effect in the nighttime LF (low frequency 30-300 kHz) signal for the Fukushima earthquake using a LF network configuration and a signal frequency of 40 kHz. They found the irregularities in the JJY-PTK (Petropavlovsk, Kamchatsky, Russia) path were detected 4-5 days before and 3 days after the Fukushima earthquake. The abnormalities were observed at Ito station in the Izu peninsula, Kamakura, Togane, and Katsuura in Chiba for the earthquake on November 14, 15, and 21, 2016. Besides VLF, Chowdhury et al. [40] also showed direct and indirect evidence of preseismic electromagnetic emissions associated with the same earthquake in Japan. They used the lithospheric emission in the range of ULF waves as observed from the Kakioka observatory in Japan for direct measurement, and the particle enhancement (electron) observed from the space-based National Oceanic and Atmospheric Administration (NOAA) satellites for indirect investigation. They observed particle bursts (PB) on the day of the Fukushima EQ. The abnormal increment in lithospheric radiation and an ionospheric depression were observed for 10 days (November 11, 2016) and 6 days (November 15, 2016) in ULF before the earthquake. The hypothesis of AGW excitation due to seismic events has been reported by Garmash et al. [12], Linkov et al. [13], and Shalimov [14]. Murayama et al. [41] utilized the middle and upper atmosphere (MU) radar to calculate the annual variation of AGW energies in Japan. Korepanov et al. [8] with the assistance of surface environmental pressing factors and attractive field information reported that AGW can be a significant boundary in the seismo-ionospheric study. Zhang et al. [42] and Yang et al. [43] utilized the sounding of the atmosphere using broadband emission radiometry (SABER) instrument onboard the thermosphere ionosphere mesosphere energetics and dynamics (TIMED) satellite temperature profile to study the potential energy (Ep) associated with AGW. Nakamura et al. [44] reported a comparative study to track down the relating seismogenic impact of some earthquakes. For the 2004, Niigata-Chuetsu earthquake (M=6.8), wavelet investigation of those boundaries, shows variances in the period of 10 to 100 minutes, which is within the scope of AGW. Yang [45] first reported that the AGW hypothesis can be used as an earthquake precursor by using the ERA5 instrument temperature profile for the 2016 Kumamoto earthquake. They observed an enhancement in AGW activity 4-6 days before the earthquake. After that, in Yang and Hayakawa [46], they reported the same hypothesis for the Tohoku earthquake (2011) using both the SABER and ERA5 temperature profiles and compared the results with the Kumamoto earthquake. In a previous study, Biswas et al. [47] reported that geomagnetic storms that occurred around the Imphal earthquake in 2016 could not contaminate the stratospheric AGW as computed from SABER. This is also a verified AGW excitation as obtained from the VLF signal. Piersanti et al. [48] detected increased AGW activity on the day of the 2018 Bayan earthquake using ERA5 data. Carbone et al. [49] computed the vertical atmospheric temperature profile using a mathematical model of the lithosphere-atmosphere interaction for seismic events and compared the findings to observations. Sasmal et al. [25] reported abnormal AGW activity 6 days before the 2020 Samos earthquake by using the SABER temperature profile. Kundu et al. [50] identified the anomalous activity in AGW around 2-20 days before the multiple earthquakes.The evidence of AGW can also be computed from the nighttime VLF/LF signal by Fourier or wavelet analysis [27, 29, 36, 51–55]. The periodic thrust produced by the gravity waves can affect the lower ionospheric characteristics by modulating the electron density profile and effective reflection height of VLF signals. This modulates the electrical conductivity in the lower ionosphere and can detect periodic changes in the amplitude and phase of VLF signals [56].Total electron content (TEC) is also found to be a potential parameter to detect seismogenic irregularities in the middle and upper ionospheres [23, 57]. It is defined as the total number of free thermal electrons in the path from a satellite to a receiver. The computation of TEC is mainly made from GPS data. During the propagation of GPS signals through the ionosphere, a group-path delay occurs due to the dispersed nature of the ionosphere. The group-path delay produces a timing error that is directly proportional to TEC [58]. The ionospheric TEC variation during an earthquake was first studied by Liu et al. [57] using a GPS receiver for the 1999 Chi-chi earthquake, having magnitude M=7.6. They found that there was a decrement in TEC during the afternoon period on days 1, 3, and 4 before this earthquake. Liu et al. [23] used a statistical method to analyses the ionospheric anomalies using the global ionosphere map (GIM) TEC during the 20 M ≥6.0 earthquakes in Taiwan from September 1999 to December 2002. After that, more similar works were studied by using GIM-TEC and GPS-TEC to study ionospheric anomalies before large earthquakes, and they found anomalies in TEC 0-17 days prior to the earthquakes [24, 59]. Using spectral and statistical analysis, Oikonomou et al. [60] observed wave-like structures having periods of 20 min and 2-5 min before the M=7.8 Nepal and M=8.3 Chile EQs in 2015. Sharma et al. [61] studied around 160 earthquakes that occurred from 2012-2018 having magnitude M≥6 using CORS (cross-origin resource sharing) installed in the north eastern states of India. They found the anomalous TEC 0-17 days prior to the earthquake and also found at least one anomalous day that is not influenced by any other solar event for each earthquake. Sasmal et al. [25] found the anomalous TEC variation around 1, 6, and 9 days prior to the earthquake using the statistical method and also found the wave-like structure, or AGW, 11 days before the earthquake obtained from GNSS-TEC using wavelet analysis.In this current manuscript, we choose the 2016 Fukushima earthquake (EQ) which took place at southeast of Namie, Fukushima Prefecture (geographic location of epicenter: 37.392° N, 141.403° E) in Japan with a magnitudeM=6.9 and a depth of 11.4 km on November 22 at 05 : 59 local time (Nov 21, 20 : 59 UTC). We use a Japanese VLF network operated by the Hayakawa Institute of Seismo Electromagnetics Co. Ltd. (Hi-SEM) and study the VLF nighttime fluctuation for the five VLF receiving locations for the transmitter JJI (22.2 kHz). These types of electromagnetic anomalies in the LAIC mechanism are corroborated by ionospheric TEC as recorded by the GNSS-IGS stations. The wavelike structures (AGW) in the VLF amplitude are studied by wavelet analysis, and the intensification of AGW is validated from the direct measurement of the SABER/TIMED outcomes. We see some new consequences in this study. We find the first evidence of altitudinal inhomogeneity in ionospheric perturbation due to electromagnetic channel. We observe almost no preseismic effect in the lower ionosphere with the help of VLF, but the previous study suggested the effect at higher altitudes with LF [39]. This study highlights the necessity of a multichannel approach to come to any definite conclusion. We see no effect in the electromagnetic channel by using VLF, but the same tool gives strong signatures of the preseismic AGW effect. It also indicates the independent working phenomena of two channels. The plan of this paper is the following. In Section 2, we explain the methodologies for the data analysis techniques. In Section 3, we present the results and discuss them. Finally, in Section 4, we conclude our findings.


## 2. Data and Methodologies

### 2.1. Analysis of VLF Radio Signal

We used multipath signal analysis around the epicenter of the earthquake in Japan. We consider a VLF network with 5 receiving stations and a transmitter JJI with a signal frequency of 22.2 kHz. The receivers are NSB (Nakashibetsu), KTU (Katsuura) in Chiba prefecture, IMZ (Imizu) in Toyama prefecture, TYH (Toyohashi) in Aichi prefecture, and ANA (Anan) in Tokushima prefecture. The locations of the VLF transmitter, receivers, the length of the great circle paths between them, and the distance from the earthquake epicenter are tabulated in Table1. We chose a ±15-day time frame from the Fukushima EQ (November 21, 2016) for our study. On November 11, 2016, another earthquake of magnitude M=6.1 occurred, with the epicenter (38.497° N, 141.566° E.), 123 km away from the Fukushima earthquake. In the VLF network, all the receiving locations use similar setups for the data acquisition where the system consists of an e-field antenna and MSK (minimum shift keying) modules. Both amplitude and phase information are recorded with a one second timing resolution. The TEC data are estimated using conventional methods as used in the GNSS-IGS system. We use MIZU (39.1351694° N, 141.1328278° E), Mizusawa, Japan, IGS station to study TEC variation. We use small-scale fluctuations in the TEC profile to check the wavelike structure at the same station. The GPS RINEX (receiver independent exchange format) observation and navigation files for the computation of TEC are taken from the IGS data archive (https://cddis.gsfc.nasa.gov/).Table 1
Locations of the VLF transmitter and receiver along with the distance from the earthquake epicenter.
TransmitterReceiverLocationDistance from the transmitter (km)Distance from the earthquake epicenter (km)JJI (22.1 kHz)32.04° N, 130.81° EANA33.89° N, 134.66° E413.81722IMZ36.78° N, 137.07° E777.84389.5KTU35.15° N, 140.31° E943.26267NSB43.54° N, 144.97° E1780.02747.8TYH34.74° N, 137.37° E677.4466The computation and calibration of TEC are done by the well-known method prescribed by Gopi Seemala by using the biasing error [62–65]. The operational region of the IGS station is approximately 300 km in radius, and it contains the epicenter of the earthquake, allowing for detecting changes in the ionosphere as well as the atmosphere. This enables us to study the condition of seismogenic impressions in the TEC profile. We have taken the geomagnetic indices such as Kp, Dst, and ap data from the World Data Center of Geomagnetism, Kyoto (http://wdc.kugi.kyoto-u.ac.jp). In Figure 1, we illustrate the locations of the VLF network and earthquake epicenters at the IGS station. The 5th Fresnel zone (FZ) [7] of all the path, along with the earthquake preparation zone (EPZ) [66] and critical zone (CZ) [67], are shown in Figure 1. Figure 1 shows that the maximum portions of all the paths are inside the EPZ, and only a small portion of the JJI-NSB path Fresnel zone is inside the CZ. The JJI-KTU path Fresnel zone is very close to the CZ. The great circle paths of the JJI to different locations are also shown in Figure 1. The IGS station MIZU is also shown within it, and it is situated inside the CZ. Figure 2 shows the geomagnetic indices during the observable period. According to Lowea and Prolss [68], the minimum value of the Dst index is employed as a criterion to classify the strengths of magnetic storms. It is to be mentioned that a moderate geomagnetic storm occurred on November 10, 2016, and it attains a minimum Dst value of −59 nT at 17 : 00. Thus, the overall situation can be treated as geomagnetically quiet, as there was no such significant contamination due to solar-terrestrial interaction. To check the anomalies in the VLF signal, we use the conventional “nighttime fluctuation” approach as described in [31]. This deals with the presence of abnormal nighttime signal amplitude shifts from the unperturbed average value during this earthquake. The trend and dispersion are computed using the following equations
(1)TrendTr.=∫t=Nbt=NeAt−AavgtdtNe−Nb,DispersionDisp.=∫t=Nbt=NeAt−Aavgt−Tr2dt,Figure 1
The red square box indicates the transmitter JJI, and straight and elliptical curves indicate the receiver locations, great circle paths, and the fifth Fresnel zones. The black circle represents the EPZ, and the red circle represents the CZ. The black diamond represents the epicenter of the earthquake, and the blue solid square indicates the IGS station MIZU.Figure 2
Dst (top panel),Kp (middle panel), and ap (bottom panel) indices activity has been shown around 15 days from the Fukushima earthquake day. ‘0’ of x-axis (blue line indicator) is the earthquake event day (November 21, 2016).Here,Nb and Ne are the start and end times of the nighttime data in seconds. Aavgt is the average amplitude at a time (t) of 31 nights, and At is the signal amplitude at a particular time. We normalized these two factors by dividing each value by the standard deviation (σ) of the total 31-night counts. According to the earlier research [6, 31], the anomalous VLF signal caused by seismic effect is characterized by a decrease in trend below the −2σ level, and an increase in dispersion over the 2σ level. To investigate the wave-like structures in the VLF nighttime signal amplitude, we use both the fast Fourier transform (FFT) and the Morlet wavelet power spectrum (WPS) analysis [69] to determine the seismogenic AGW. We performed this for 15 days, from November 6 to 21, 2016. For the FFT, we use the conventional rectangular window. To compute the wavelet power spectrum (WPS), we first rearrange the data in a one-minute time sample and then subtract each amplitude value from its ten-minute running mean, which gives the fluctuation in amplitude. We draw the cone of influence (CoI) in WPS. The values beyond this CoI zone are not deemed legitimate due to the addition of zeros (zero padding) required to convert the total number of data points to a power of two to complete the WPS.


### 2.2. AGW Computation Using SABER/TIMED

SABER is an onboard instrument of the TIMED satellite. The satellite was launched on December 7, 2001, into an orbit with a height of 625 km. The period of the satellite is 1.7 hours with an inclination at 74.1° [70]. SABER made its first observation in January 2002. It obtains the temperature with an altitude coverage of 100 km by using wavelengths ranging from 1.27 to 17 μm. From northward viewing, the latitudinal range of SABER extends between 50° S and 82° N, and from southward view it is 50° N and 82° S. The observation viewing technique of SABER changes every 60 days. Remsberg et al. [71] discussed the method of retribution of temperature from SABER. The mechanism to obtain the AGW from the temperature profile has been obtained by several methods in the past few decades [25, 46, 48, 72–76]. We used the following mechanism to extract the AGW from the SABER temperature profile. At first, we take the temperature profile for the altitude range 20–100 km for the region of our study from the SABER archive data http://saber.gats-inc.com/, and we take the logarithm of the obtained individual temperature profiles. Then, a third-order polynomial is fitted on the logarithm temperature profile, and we extract the fitted value. To get the residuals, we subtract the fitted profile from the original one. As AGW has a wavelength greater than 4 km, to remove the other waves having small wavelengths, a 4 km boxcar filter is applied to the residuals of the individual profiles. After filtration, the filtered data is added back to the fitted profile, and we got the final profile. An antilogarithm of the final profiles is known as the least square fit (LSF) which are used to obtain the daily zonal mean temperature and other zonal wave components from 1–5. The background temperature (T0) is obtained from the summation of all wave components 0–5. The background temperature profiles are subtracted from the original one to obtain the perturbation temperature (Tp). The obtained background temperature profiles corresponding to the altitude are put in Equation (2) to get the Brunt Väisälä frequency (N) of the respective profiles as written as,
(2)N2=gT0∂T0∂z+gcp,where z is the altitude, cp is the specific heat at constant pressure, and g is the acceleration due to gravity.The potential energy (Ep) associated with the AGW for individual temperature profile can be obtained by putting the values of perturbation and background temperature in Equation (3) as,
(3)Ep=12gN2TpT02,N is the Brunt Vaisala frequency. This method is already used in our previous works to study the activity of AGW during the 2016 Imphal earthquake, [47] and the 2020 Samos earthquake [25]. The method is originally used by [74] to extract the GW from the SOFIE temperature profile to study the activity of AGW during multiple sudden stratospheric warming (SSW) events. To study the AGW activity during the EQ, we choose 29 days, from November 1 to 29, 2016. The region for this study has a span between 25°-50° N (latitude range) and 125°-150° E (longitude range).


### 2.3. Ionospheric Perturbation from GPS-TEC

We use GPS-TEC data to obtain the AGW which is nothing but a traveling ionospheric disturbance (TID). The disturbances created by the TIDs are small-scale fluctuations that can be generated due to earthquakes and meteorological phenomena. This small-scale fluctuation can be captured by a GPS signal during its propagation. To detect the anomaly, we compute the medianX by using the past 15 days’ vertical TEC (VTEC) values, the associated interquartile range IQR, and the upper bound (UB), and lower bound (LB) at a certain UTC. We use the same method for the computation of anomaly in VTEC that is used by Liu et al. [23]. Under the assumption of a normal distribution with standard deviation (σ) for the VTEC, the expected value of IQR is 1.34σ [75]. After computation of the upper bound and lower band, we calculate the anomaly. The anomaly is defined as the enhancement of diurnal TEC concerning the upper bound or decrement of diurnal TEC for the lower bound with the confidence level of 80–85% [24].


### 2.4. AGW Computed from GPS-TEC

We use a method to get small-scale fluctuation from the TEC data. We use the MATLAB algorithm “Savitzky-Golay filtering” (sgolayfilt)’ to filter or smooth the TEC values. We use fifth-order sgolayfilt with a frame length or time window of 90 minutes. We repeat the process two times to remove the other noises from the data, and we get the fluctuation. Fitted data is subtracted from the observed data to get small-scale fluctuations (Equation (4)) [25, 76]. If the small-scale fluctuation is dVTEC, the VTEC is the observed TEC profile, and the VTECf is the fitted profile, then it can be written as,
(4)dVTEC=VTEC−VTECf.The fitting technique “Savitzky-Golay filtering” (sgolayfilt) [77, 78] process is expressed in Equation (5). Savitzky and Golay [77] demonstrated a collection of integers (B−n, B−(n1)..., Bn−1, Bn), which may be used to calculate the weighting coefficients for a smoothing operation. The use of these weighting coefficients, also known as convolution integers, is identical to fitting the data to a polynomial, as stated, but is more computationally efficient and quicker. The Savitzky-Golay method produces the smoothed data point (zk)s as,
(5)zks=∑j=−nj=nBkzk+1∑k=−nk=nBk.After getting a small-scale fluctuations ordVTEC, the spectral analysis of dVTECis obtained to get the wavelike structures. We do a wavelet analysis of dVTEC by using the complex Morlet continuous wavelet analysis technique [69]. The wavelike structure generated within the CoI can be treated as seismogenic AGW [25, 76].


## 2.1. Analysis of VLF Radio Signal

We used multipath signal analysis around the epicenter of the earthquake in Japan. We consider a VLF network with 5 receiving stations and a transmitter JJI with a signal frequency of 22.2 kHz. The receivers are NSB (Nakashibetsu), KTU (Katsuura) in Chiba prefecture, IMZ (Imizu) in Toyama prefecture, TYH (Toyohashi) in Aichi prefecture, and ANA (Anan) in Tokushima prefecture. The locations of the VLF transmitter, receivers, the length of the great circle paths between them, and the distance from the earthquake epicenter are tabulated in Table1. We chose a ±15-day time frame from the Fukushima EQ (November 21, 2016) for our study. On November 11, 2016, another earthquake of magnitude M=6.1 occurred, with the epicenter (38.497° N, 141.566° E.), 123 km away from the Fukushima earthquake. In the VLF network, all the receiving locations use similar setups for the data acquisition where the system consists of an e-field antenna and MSK (minimum shift keying) modules. Both amplitude and phase information are recorded with a one second timing resolution. The TEC data are estimated using conventional methods as used in the GNSS-IGS system. We use MIZU (39.1351694° N, 141.1328278° E), Mizusawa, Japan, IGS station to study TEC variation. We use small-scale fluctuations in the TEC profile to check the wavelike structure at the same station. The GPS RINEX (receiver independent exchange format) observation and navigation files for the computation of TEC are taken from the IGS data archive (https://cddis.gsfc.nasa.gov/).Table 1
Locations of the VLF transmitter and receiver along with the distance from the earthquake epicenter.
TransmitterReceiverLocationDistance from the transmitter (km)Distance from the earthquake epicenter (km)JJI (22.1 kHz)32.04° N, 130.81° EANA33.89° N, 134.66° E413.81722IMZ36.78° N, 137.07° E777.84389.5KTU35.15° N, 140.31° E943.26267NSB43.54° N, 144.97° E1780.02747.8TYH34.74° N, 137.37° E677.4466The computation and calibration of TEC are done by the well-known method prescribed by Gopi Seemala by using the biasing error [62–65]. The operational region of the IGS station is approximately 300 km in radius, and it contains the epicenter of the earthquake, allowing for detecting changes in the ionosphere as well as the atmosphere. This enables us to study the condition of seismogenic impressions in the TEC profile. We have taken the geomagnetic indices such as Kp, Dst, and ap data from the World Data Center of Geomagnetism, Kyoto (http://wdc.kugi.kyoto-u.ac.jp). In Figure 1, we illustrate the locations of the VLF network and earthquake epicenters at the IGS station. The 5th Fresnel zone (FZ) [7] of all the path, along with the earthquake preparation zone (EPZ) [66] and critical zone (CZ) [67], are shown in Figure 1. Figure 1 shows that the maximum portions of all the paths are inside the EPZ, and only a small portion of the JJI-NSB path Fresnel zone is inside the CZ. The JJI-KTU path Fresnel zone is very close to the CZ. The great circle paths of the JJI to different locations are also shown in Figure 1. The IGS station MIZU is also shown within it, and it is situated inside the CZ. Figure 2 shows the geomagnetic indices during the observable period. According to Lowea and Prolss [68], the minimum value of the Dst index is employed as a criterion to classify the strengths of magnetic storms. It is to be mentioned that a moderate geomagnetic storm occurred on November 10, 2016, and it attains a minimum Dst value of −59 nT at 17 : 00. Thus, the overall situation can be treated as geomagnetically quiet, as there was no such significant contamination due to solar-terrestrial interaction. To check the anomalies in the VLF signal, we use the conventional “nighttime fluctuation” approach as described in [31]. This deals with the presence of abnormal nighttime signal amplitude shifts from the unperturbed average value during this earthquake. The trend and dispersion are computed using the following equations
(1)TrendTr.=∫t=Nbt=NeAt−AavgtdtNe−Nb,DispersionDisp.=∫t=Nbt=NeAt−Aavgt−Tr2dt,Figure 1
The red square box indicates the transmitter JJI, and straight and elliptical curves indicate the receiver locations, great circle paths, and the fifth Fresnel zones. The black circle represents the EPZ, and the red circle represents the CZ. The black diamond represents the epicenter of the earthquake, and the blue solid square indicates the IGS station MIZU.Figure 2
Dst (top panel),Kp (middle panel), and ap (bottom panel) indices activity has been shown around 15 days from the Fukushima earthquake day. ‘0’ of x-axis (blue line indicator) is the earthquake event day (November 21, 2016).Here,Nb and Ne are the start and end times of the nighttime data in seconds. Aavgt is the average amplitude at a time (t) of 31 nights, and At is the signal amplitude at a particular time. We normalized these two factors by dividing each value by the standard deviation (σ) of the total 31-night counts. According to the earlier research [6, 31], the anomalous VLF signal caused by seismic effect is characterized by a decrease in trend below the −2σ level, and an increase in dispersion over the 2σ level. To investigate the wave-like structures in the VLF nighttime signal amplitude, we use both the fast Fourier transform (FFT) and the Morlet wavelet power spectrum (WPS) analysis [69] to determine the seismogenic AGW. We performed this for 15 days, from November 6 to 21, 2016. For the FFT, we use the conventional rectangular window. To compute the wavelet power spectrum (WPS), we first rearrange the data in a one-minute time sample and then subtract each amplitude value from its ten-minute running mean, which gives the fluctuation in amplitude. We draw the cone of influence (CoI) in WPS. The values beyond this CoI zone are not deemed legitimate due to the addition of zeros (zero padding) required to convert the total number of data points to a power of two to complete the WPS.


## 2.2. AGW Computation Using SABER/TIMED

SABER is an onboard instrument of the TIMED satellite. The satellite was launched on December 7, 2001, into an orbit with a height of 625 km. The period of the satellite is 1.7 hours with an inclination at 74.1° [70]. SABER made its first observation in January 2002. It obtains the temperature with an altitude coverage of 100 km by using wavelengths ranging from 1.27 to 17 μm. From northward viewing, the latitudinal range of SABER extends between 50° S and 82° N, and from southward view it is 50° N and 82° S. The observation viewing technique of SABER changes every 60 days. Remsberg et al. [71] discussed the method of retribution of temperature from SABER. The mechanism to obtain the AGW from the temperature profile has been obtained by several methods in the past few decades [25, 46, 48, 72–76]. We used the following mechanism to extract the AGW from the SABER temperature profile. At first, we take the temperature profile for the altitude range 20–100 km for the region of our study from the SABER archive data http://saber.gats-inc.com/, and we take the logarithm of the obtained individual temperature profiles. Then, a third-order polynomial is fitted on the logarithm temperature profile, and we extract the fitted value. To get the residuals, we subtract the fitted profile from the original one. As AGW has a wavelength greater than 4 km, to remove the other waves having small wavelengths, a 4 km boxcar filter is applied to the residuals of the individual profiles. After filtration, the filtered data is added back to the fitted profile, and we got the final profile. An antilogarithm of the final profiles is known as the least square fit (LSF) which are used to obtain the daily zonal mean temperature and other zonal wave components from 1–5. The background temperature (T0) is obtained from the summation of all wave components 0–5. The background temperature profiles are subtracted from the original one to obtain the perturbation temperature (Tp). The obtained background temperature profiles corresponding to the altitude are put in Equation (2) to get the Brunt Väisälä frequency (N) of the respective profiles as written as,
(2)N2=gT0∂T0∂z+gcp,where z is the altitude, cp is the specific heat at constant pressure, and g is the acceleration due to gravity.The potential energy (Ep) associated with the AGW for individual temperature profile can be obtained by putting the values of perturbation and background temperature in Equation (3) as,
(3)Ep=12gN2TpT02,N is the Brunt Vaisala frequency. This method is already used in our previous works to study the activity of AGW during the 2016 Imphal earthquake, [47] and the 2020 Samos earthquake [25]. The method is originally used by [74] to extract the GW from the SOFIE temperature profile to study the activity of AGW during multiple sudden stratospheric warming (SSW) events. To study the AGW activity during the EQ, we choose 29 days, from November 1 to 29, 2016. The region for this study has a span between 25°-50° N (latitude range) and 125°-150° E (longitude range).


## 2.3. Ionospheric Perturbation from GPS-TEC

We use GPS-TEC data to obtain the AGW which is nothing but a traveling ionospheric disturbance (TID). The disturbances created by the TIDs are small-scale fluctuations that can be generated due to earthquakes and meteorological phenomena. This small-scale fluctuation can be captured by a GPS signal during its propagation. To detect the anomaly, we compute the medianX by using the past 15 days’ vertical TEC (VTEC) values, the associated interquartile range IQR, and the upper bound (UB), and lower bound (LB) at a certain UTC. We use the same method for the computation of anomaly in VTEC that is used by Liu et al. [23]. Under the assumption of a normal distribution with standard deviation (σ) for the VTEC, the expected value of IQR is 1.34σ [75]. After computation of the upper bound and lower band, we calculate the anomaly. The anomaly is defined as the enhancement of diurnal TEC concerning the upper bound or decrement of diurnal TEC for the lower bound with the confidence level of 80–85% [24].


## 2.4. AGW Computed from GPS-TEC

We use a method to get small-scale fluctuation from the TEC data. We use the MATLAB algorithm “Savitzky-Golay filtering” (sgolayfilt)’ to filter or smooth the TEC values. We use fifth-order sgolayfilt with a frame length or time window of 90 minutes. We repeat the process two times to remove the other noises from the data, and we get the fluctuation. Fitted data is subtracted from the observed data to get small-scale fluctuations (Equation (4)) [25, 76]. If the small-scale fluctuation is dVTEC, the VTEC is the observed TEC profile, and the VTECf is the fitted profile, then it can be written as,
(4)dVTEC=VTEC−VTECf.The fitting technique “Savitzky-Golay filtering” (sgolayfilt) [77, 78] process is expressed in Equation (5). Savitzky and Golay [77] demonstrated a collection of integers (B−n, B−(n1)..., Bn−1, Bn), which may be used to calculate the weighting coefficients for a smoothing operation. The use of these weighting coefficients, also known as convolution integers, is identical to fitting the data to a polynomial, as stated, but is more computationally efficient and quicker. The Savitzky-Golay method produces the smoothed data point (zk)s as,
(5)zks=∑j=−nj=nBkzk+1∑k=−nk=nBk.After getting a small-scale fluctuations ordVTEC, the spectral analysis of dVTECis obtained to get the wavelike structures. We do a wavelet analysis of dVTEC by using the complex Morlet continuous wavelet analysis technique [69]. The wavelike structure generated within the CoI can be treated as seismogenic AGW [25, 76].


## 3. Results

### 3.1. Results Obtained from VLF Signal

The normalized trend and dispersion as computed from the VLF nighttime fluctuations at 5 different receiving stations (ANA, IMZ, KTU, NSB, and TYH) are shown in Figures3 and 4. All the anomalous days are marked with red color. Apart from all the stations, we notice a significant reduction in the trend only for the IMZ station below the −2σ threshold on 6 to 7 days before the EQ (Figure 3(b)). No such significant depletion in the trend profile is observed for the rest of the stations. In the postearthquake period, the trend decreases below 2σ, after 2 days for the ANA station. As it is evident from Figure 4 for the dispersion, a fluctuation of more than 2σ is detected in the NSB station one day before the earthquake. For the TYH station, there are two enhancements in the dispersion value. The initial peak is observed two weeks before the earthquake, and that is four days before the first earthquake on November 11. Postevent fluctuation observed only in ANA above 4σ level, after 6 days of the earthquake. This can be associated with the numerous aftershocks.Figure 3
Normalized trend obtained from (a) ANA, (b) IMZ, (c) KTU, (d) NSB, and (e) TYH stations. ‘0’ in thex-axis is the date of the earthquake (November 21, 2016) and the dotted horizontal lines indicates the σ levels.
(a)(b)(c)(d)(e)Figure 4
Same as Figure3, but for normalized dispersion.
(a)
ANA(b)
IMZ(c)
KTU(d)
NSB(e)
TYHBy using the FFT and wavelet analysis, we demonstrate the presence of wavelike structures in the nighttime data in Figures5 and 6. Figure 5 shows the normalized FFT periodic spectrums of (a) ANA, (b) IMZ, (c) KTU, (d) NSB, and (e) TYH stations. The red lines denote the days with waves having periodicity similar to AGW have enhanced intensity. For the five stations, we find (a) November 7, 8, and 9; (b) November 12, 13, 15, and 16; (c) November 12; (d) November 11 and 21; and (e) November 7, 12, 17, and 18 can be treated as seismically active days when there is the presence of AGW intensification. The black curves denote the rest of the other quiet days.Figure 5
Normalized FFT periodic spectrums of (a) ANA, (b) IMZ, (c) KTU, (d) NSB, and (e) TYH stations. Horizontal axes refer to periodicity in minutes up to 128 minutes and vertical axes refer to the Fourier amplitude in dB. Red lines are denoting the higher periodicity activity days. The wavelike structures are found to be having the periods ranging from 40 to 100 minutes.
(a)(b)(c)(d)(e)Figure 6
Wavelet power spectrum of (a) ANA, (b) IMZ, (c) KTU, (d) NSB, and (e) TYH stations. Horizontal axis refers to time of the night after 12 : 00 : 00 UTC. The white parabolic line is the cone of influence (CoI). Vertical axis refers to periodicity in minutes and power is the magnitude of periodic structure.
(a)(b)(c)(d)(e)The outcomes of FFT analysis give a mixed impression. For ANA, (Figure5(a)), we detect larger amplitude of wave-like structures of periodicity ~82 minutes on November 7, ~85 minutes on November 8, and 70 minutes on November 9, 2016. From the IMZ station, waves of intermediate amplitude are observed having periodicity ~55–60 minutes (Figure 5(b)) on November 12, 13, 14, and 16. On November 15, there are two peaks observed ~60 minutes and ~106 minutes but the second peak is not so sharp. In the case of KTU, the periodic structure presents only on November 12 (Figure 5(c)) with a periodicity ~39, 58, and 83 minutes. From NSB, the periodic structures are not so prominent (Figure 5(d)). It is difficult to determine the prominent periodicity on November 11 and 21 due to the lack of sharp peaks. For TYH (Figure 5(e)), we observe waves of having periods ~105 minutes on November 7, ~49 and 69 minutes on November 12, ~56 minutes on November 17, and~58 minutes on November 18, 2016.The wavelet power spectrum of the ANA station (Figure6(a)) indicates that an intense wave-like structure appeared within the CoI on November 8 with a period of 65-90 minutes. Another intense wave-like structure is observed on November 18 with a period of ∼48 minutes. For IMZ (Figure 6(b)), the WPS outcomes match with its FFT results. Periodicity of around 64 minutes is obtained on November 12, 13, 15, and 16. In Figure 6(c), one can see that for KTU, intensification in AGW is observed with a period of around 40 minutes on November 12. Waves of periods ~35-50 minutes are observed for NSB station (Figure 6(d)) on November 11 and 21. From the data of TYH (Figure 6(e)), intense AGW activities are observed on November 14 with a period around ∼35 minutes, on November 17 with periodicity around 40-60 minutes, and on November 18 with the period around ~35–64 minutes, and also less intense periodic structures have been found on November 07 (periodicity around ~105 minutes, also found in FFT, and periodicity around ~32 minutes).


### 3.2. Observed Results from SABER

After the computation of Ep, a nine-dimensional matrix is obtained for individual temperature profiles. The nine-dimensional matrix consists of latitude, longitude, date or day of the year in UT, altitude, original SABER temperature profile, reconstructed fitted temperature profile, perturbation temperature, Brunt Vaisala frequency, and AGW-associated potential energy (Ep). At first, we obtain the altitude profiles of the above-said parameters. The significant enhancement in Ep associated with AGW activity is observed from November 16 to 19, 2016. The maximum activity of AGW is observed on November 19, 2016. Figure7 shows that a maximum value of Ep on this day is 11.47 J/kg at an altitude of 44 km.Figure 7
Altitude variation of (a) temperature (black solid line) in Kelvin with 3rd order polynomial fitted value in red dashed line, (b) perturbation temperature T J in solid black line and the nonperturbed temperature indicated in blue dashed line, (c) Brunt Vaisala frequency N2 in rad/S2, and (d) calculated potential energy Ep in J/kg from 30 to 50 km on November 19, 2016 for Fukushima EQ.The temporal variation or time-altitude variation of Ep with the associated AGW activity is also obtained from the nine-dimensional matrix. From Figure8, an enhancement of Ep is observed from November 16 to 19, 2016 around 44 to 46 km altitude. The AGW activity is significantly enhanced just 2-4 days before the EQ on November 21, 2016.Figure 8
Variation of Ep during November 1 to 29, 2016 for Fukushima EQ. Along X and Y axis, we present the date and altitude in km, respectively. The black dashed line indicates the earthquake day.The spatial variation of Ep is also obtained to observe the AGW activity or enhancement of Ep near the epicenter of the earthquake. The spatial variation is obtained around 46 km altitude during maximum activity of AGW time which is already obtained from the altitude and time-altitude variation. From Figure9, it is observed that for 46 km there is a high Ep enhancement around the epicenter during November 16–19, 2016. Maximum enhancement is observed on the 18 of November around the epicenter. The enhancement in Ep associated with seismogenic AGW is found to be prominent around 44 to 46 km altitude 2 to 4 days (November 16 to 19, 2016) prior to the earthquake. The maximum enhancement is observed on November 18, 2016, 528 at 46 km altitude around the epicenter. The enhancement of Ep is started on November 16 and traveled from the west to east side of the epicenter. We also observe enhancement on 24, 26, and 28 November around the epicenter that is possibly connected with the numerous aftershocks as reported by USGS.Figure 9
Spatial variation of potential energy at an altitude of 46 km from November 16 to November 30, 2016 for the 2016 Fukushima earthquake. Along thex-axis and y-axis, we present the longitude and latitude, respectively. The magenta diamond represents the epicenter of the earthquake and the white lines indicate the country border.


### 3.3. Results Obtained from GPS-TEC

The diurnal variation of VTEC from November 6 to November 23, 2016, along with the upper bound and lower are shown in the upper panel of Figure10. The enhancement in VTEC is much higher than the upper bound on November 11 and 15, 2016. In these two days, the enhancement in TEC is 4.5 and 2.9 TECU, respectively. The change in TEC is shown in the lower panel of Figure 10. It is evident that the anomaly is observed in VTEC around 6 and 10 days before the earthquake. But November 11 is the second day after beginning the moderate magnetic storm on November 10 (Dstmin=−59 nT, November 10, 17 UT). Thus, on the November 11, TEC anomaly is most likely associated with this moderate magnetic storm. As November 15 is geomagnetically quiet; so the anomaly observed on November 15, 2016, can be associated with the earthquake.Figure 10
(a) Diurnal variation of VTEC from November 6 to November 23, 2016 for the Fukushima earthquake along with the upper bound and lower bound. The black line indicates the TEC variation and the red and green line indicate the upper and lower bound, respectively. (b) Change in VTEC for the same time interval.
(a)(b)To verify the AGW excitation as obtained from SABER, wavelike structures are computed from the small-scale fluctuations in the VTEC outcomes. The normal unperturbed ionospheric condition is found to be−0.3≥dVTEC≤0.3 [25, 76]. To assign a small-scale fluctuation as a seismogenic perturbation, the value of dVTEC must be outside the above-mentioned range. The wavelet spectrum as computed from the small-scale fluctuation satisfying the said condition is shown in Figure 11.Figure 11
Wavelet spectrum as computed from small-scale fluctuations in VTEC from November 9 to 23, 2016. The white line represents the cone of influence (CoI).From the wavelet spectrum, wave-like structures of period 60-90 minutes are observed on November 18, 2016, inside the CoI. Similar formation of AGW is also witnessed on November 9 and 21 having a period of 60-80 minutes and 60-70 minutes, respectively. The VTEC variation along with the fitted VTEC, the small-scale fluctuation, and the wavelet spectrum on November 18, 2016, is presented in Figure12. A comparison of the maximum changes in the context of the temporal weightage of all the methods is tabulated in Table 2.Figure 12
Diurnal variation of observed VTEC along with fitted curve. The black solid line represents the observed VTEC, and red dotted line represents the fitted VTEC on November 18, 2016 (upper panel). (Middle panel) Variation of dVTEC with the upper and lower normal values (red lines). The wavelet power spectrum (lower panel) with cone of influence (CoI) (white line).Table 2
Comparative study of preearthquake effect.
VLF-NTFVTECAGW from VLFdVTECAGW from VTECMaximumEpANA-no effectNov 11ANA-Nov 07-09, 18-19Nov 09Nov 09Nov 16IMZ-Nov 14-15 (trend)Nov 15IMZ-Nov 9, 11, 13, 15-16Nov 16Nov 18Nov 17KTU-no effectKTU-Nov 12Nov 18Nov 21Nov 18NSB-Nov 20 (Disp.)NSB-Nov 11, 21Nov 19Nov 19TYH-Nov 07 (Disp.)TYH-Nov 7, 14, 17-18Nov 21


## 3.1. Results Obtained from VLF Signal

The normalized trend and dispersion as computed from the VLF nighttime fluctuations at 5 different receiving stations (ANA, IMZ, KTU, NSB, and TYH) are shown in Figures3 and 4. All the anomalous days are marked with red color. Apart from all the stations, we notice a significant reduction in the trend only for the IMZ station below the −2σ threshold on 6 to 7 days before the EQ (Figure 3(b)). No such significant depletion in the trend profile is observed for the rest of the stations. In the postearthquake period, the trend decreases below 2σ, after 2 days for the ANA station. As it is evident from Figure 4 for the dispersion, a fluctuation of more than 2σ is detected in the NSB station one day before the earthquake. For the TYH station, there are two enhancements in the dispersion value. The initial peak is observed two weeks before the earthquake, and that is four days before the first earthquake on November 11. Postevent fluctuation observed only in ANA above 4σ level, after 6 days of the earthquake. This can be associated with the numerous aftershocks.Figure 3
Normalized trend obtained from (a) ANA, (b) IMZ, (c) KTU, (d) NSB, and (e) TYH stations. ‘0’ in thex-axis is the date of the earthquake (November 21, 2016) and the dotted horizontal lines indicates the σ levels.
(a)(b)(c)(d)(e)Figure 4
Same as Figure3, but for normalized dispersion.
(a)
ANA(b)
IMZ(c)
KTU(d)
NSB(e)
TYHBy using the FFT and wavelet analysis, we demonstrate the presence of wavelike structures in the nighttime data in Figures5 and 6. Figure 5 shows the normalized FFT periodic spectrums of (a) ANA, (b) IMZ, (c) KTU, (d) NSB, and (e) TYH stations. The red lines denote the days with waves having periodicity similar to AGW have enhanced intensity. For the five stations, we find (a) November 7, 8, and 9; (b) November 12, 13, 15, and 16; (c) November 12; (d) November 11 and 21; and (e) November 7, 12, 17, and 18 can be treated as seismically active days when there is the presence of AGW intensification. The black curves denote the rest of the other quiet days.Figure 5
Normalized FFT periodic spectrums of (a) ANA, (b) IMZ, (c) KTU, (d) NSB, and (e) TYH stations. Horizontal axes refer to periodicity in minutes up to 128 minutes and vertical axes refer to the Fourier amplitude in dB. Red lines are denoting the higher periodicity activity days. The wavelike structures are found to be having the periods ranging from 40 to 100 minutes.
(a)(b)(c)(d)(e)Figure 6
Wavelet power spectrum of (a) ANA, (b) IMZ, (c) KTU, (d) NSB, and (e) TYH stations. Horizontal axis refers to time of the night after 12 : 00 : 00 UTC. The white parabolic line is the cone of influence (CoI). Vertical axis refers to periodicity in minutes and power is the magnitude of periodic structure.
(a)(b)(c)(d)(e)The outcomes of FFT analysis give a mixed impression. For ANA, (Figure5(a)), we detect larger amplitude of wave-like structures of periodicity ~82 minutes on November 7, ~85 minutes on November 8, and 70 minutes on November 9, 2016. From the IMZ station, waves of intermediate amplitude are observed having periodicity ~55–60 minutes (Figure 5(b)) on November 12, 13, 14, and 16. On November 15, there are two peaks observed ~60 minutes and ~106 minutes but the second peak is not so sharp. In the case of KTU, the periodic structure presents only on November 12 (Figure 5(c)) with a periodicity ~39, 58, and 83 minutes. From NSB, the periodic structures are not so prominent (Figure 5(d)). It is difficult to determine the prominent periodicity on November 11 and 21 due to the lack of sharp peaks. For TYH (Figure 5(e)), we observe waves of having periods ~105 minutes on November 7, ~49 and 69 minutes on November 12, ~56 minutes on November 17, and~58 minutes on November 18, 2016.The wavelet power spectrum of the ANA station (Figure6(a)) indicates that an intense wave-like structure appeared within the CoI on November 8 with a period of 65-90 minutes. Another intense wave-like structure is observed on November 18 with a period of ∼48 minutes. For IMZ (Figure 6(b)), the WPS outcomes match with its FFT results. Periodicity of around 64 minutes is obtained on November 12, 13, 15, and 16. In Figure 6(c), one can see that for KTU, intensification in AGW is observed with a period of around 40 minutes on November 12. Waves of periods ~35-50 minutes are observed for NSB station (Figure 6(d)) on November 11 and 21. From the data of TYH (Figure 6(e)), intense AGW activities are observed on November 14 with a period around ∼35 minutes, on November 17 with periodicity around 40-60 minutes, and on November 18 with the period around ~35–64 minutes, and also less intense periodic structures have been found on November 07 (periodicity around ~105 minutes, also found in FFT, and periodicity around ~32 minutes).


## 3.2. Observed Results from SABER

After the computation of Ep, a nine-dimensional matrix is obtained for individual temperature profiles. The nine-dimensional matrix consists of latitude, longitude, date or day of the year in UT, altitude, original SABER temperature profile, reconstructed fitted temperature profile, perturbation temperature, Brunt Vaisala frequency, and AGW-associated potential energy (Ep). At first, we obtain the altitude profiles of the above-said parameters. The significant enhancement in Ep associated with AGW activity is observed from November 16 to 19, 2016. The maximum activity of AGW is observed on November 19, 2016. Figure7 shows that a maximum value of Ep on this day is 11.47 J/kg at an altitude of 44 km.Figure 7
Altitude variation of (a) temperature (black solid line) in Kelvin with 3rd order polynomial fitted value in red dashed line, (b) perturbation temperature T J in solid black line and the nonperturbed temperature indicated in blue dashed line, (c) Brunt Vaisala frequency N2 in rad/S2, and (d) calculated potential energy Ep in J/kg from 30 to 50 km on November 19, 2016 for Fukushima EQ.The temporal variation or time-altitude variation of Ep with the associated AGW activity is also obtained from the nine-dimensional matrix. From Figure8, an enhancement of Ep is observed from November 16 to 19, 2016 around 44 to 46 km altitude. The AGW activity is significantly enhanced just 2-4 days before the EQ on November 21, 2016.Figure 8
Variation of Ep during November 1 to 29, 2016 for Fukushima EQ. Along X and Y axis, we present the date and altitude in km, respectively. The black dashed line indicates the earthquake day.The spatial variation of Ep is also obtained to observe the AGW activity or enhancement of Ep near the epicenter of the earthquake. The spatial variation is obtained around 46 km altitude during maximum activity of AGW time which is already obtained from the altitude and time-altitude variation. From Figure9, it is observed that for 46 km there is a high Ep enhancement around the epicenter during November 16–19, 2016. Maximum enhancement is observed on the 18 of November around the epicenter. The enhancement in Ep associated with seismogenic AGW is found to be prominent around 44 to 46 km altitude 2 to 4 days (November 16 to 19, 2016) prior to the earthquake. The maximum enhancement is observed on November 18, 2016, 528 at 46 km altitude around the epicenter. The enhancement of Ep is started on November 16 and traveled from the west to east side of the epicenter. We also observe enhancement on 24, 26, and 28 November around the epicenter that is possibly connected with the numerous aftershocks as reported by USGS.Figure 9
Spatial variation of potential energy at an altitude of 46 km from November 16 to November 30, 2016 for the 2016 Fukushima earthquake. Along thex-axis and y-axis, we present the longitude and latitude, respectively. The magenta diamond represents the epicenter of the earthquake and the white lines indicate the country border.


## 3.3. Results Obtained from GPS-TEC

The diurnal variation of VTEC from November 6 to November 23, 2016, along with the upper bound and lower are shown in the upper panel of Figure10. The enhancement in VTEC is much higher than the upper bound on November 11 and 15, 2016. In these two days, the enhancement in TEC is 4.5 and 2.9 TECU, respectively. The change in TEC is shown in the lower panel of Figure 10. It is evident that the anomaly is observed in VTEC around 6 and 10 days before the earthquake. But November 11 is the second day after beginning the moderate magnetic storm on November 10 (Dstmin=−59 nT, November 10, 17 UT). Thus, on the November 11, TEC anomaly is most likely associated with this moderate magnetic storm. As November 15 is geomagnetically quiet; so the anomaly observed on November 15, 2016, can be associated with the earthquake.Figure 10
(a) Diurnal variation of VTEC from November 6 to November 23, 2016 for the Fukushima earthquake along with the upper bound and lower bound. The black line indicates the TEC variation and the red and green line indicate the upper and lower bound, respectively. (b) Change in VTEC for the same time interval.
(a)(b)To verify the AGW excitation as obtained from SABER, wavelike structures are computed from the small-scale fluctuations in the VTEC outcomes. The normal unperturbed ionospheric condition is found to be−0.3≥dVTEC≤0.3 [25, 76]. To assign a small-scale fluctuation as a seismogenic perturbation, the value of dVTEC must be outside the above-mentioned range. The wavelet spectrum as computed from the small-scale fluctuation satisfying the said condition is shown in Figure 11.Figure 11
Wavelet spectrum as computed from small-scale fluctuations in VTEC from November 9 to 23, 2016. The white line represents the cone of influence (CoI).From the wavelet spectrum, wave-like structures of period 60-90 minutes are observed on November 18, 2016, inside the CoI. Similar formation of AGW is also witnessed on November 9 and 21 having a period of 60-80 minutes and 60-70 minutes, respectively. The VTEC variation along with the fitted VTEC, the small-scale fluctuation, and the wavelet spectrum on November 18, 2016, is presented in Figure12. A comparison of the maximum changes in the context of the temporal weightage of all the methods is tabulated in Table 2.Figure 12
Diurnal variation of observed VTEC along with fitted curve. The black solid line represents the observed VTEC, and red dotted line represents the fitted VTEC on November 18, 2016 (upper panel). (Middle panel) Variation of dVTEC with the upper and lower normal values (red lines). The wavelet power spectrum (lower panel) with cone of influence (CoI) (white line).Table 2
Comparative study of preearthquake effect.
VLF-NTFVTECAGW from VLFdVTECAGW from VTECMaximumEpANA-no effectNov 11ANA-Nov 07-09, 18-19Nov 09Nov 09Nov 16IMZ-Nov 14-15 (trend)Nov 15IMZ-Nov 9, 11, 13, 15-16Nov 16Nov 18Nov 17KTU-no effectKTU-Nov 12Nov 18Nov 21Nov 18NSB-Nov 20 (Disp.)NSB-Nov 11, 21Nov 19Nov 19TYH-Nov 07 (Disp.)TYH-Nov 7, 14, 17-18Nov 21


## 4. Discussion

The simultaneous study of ground and space-based methods using two major channels of the LAIC mechanism gives overall satisfactory outcomes. However, we observed predominant effects in the acoustic channel (AGW) in comparison to the electromagnetic channel (TEC and VLF waves). We observe relatively less prominent seismogenic effects in the VLF nighttime fluctuations. Out of five stations, we detect a significant anomalous trend in the IMZ signal 6 to 7 days before the earthquake, as well as anomalous dispersion in the NSB signal on 1 day and 14 days before the earthquake. It is also evident that for NSB and TYH stations, preseismic depletion in trend is observed within the 10 days before the earthquake, but the variations are not so prominent (within 1.5 to 2σ levels). The ANA and NSB stations are situated far from the earthquake epicenter. Despite the fact that the distance from JJI to ANA is shorter in comparison to that for the NSB station, the ANA signal does not present any preseismic effect in the trend and dispersion. In spite of having a larger distance from JJI, the NSB signal shows a significant increase in dispersion with a moderate depletion in trend. The JJI-NSB path has clear overlapping (Figure 1) with the EPZ [66] which could be the possible reason behind it.The IMZ and KTU stations having intermediate path lengths (Table1) show a mixed nature in anomalies. IMZ shows a strong depletion in trend and a medium enhancement in dispersion, whereas KTU shows only a medium enchantment in dispersion. Sitting a bit closer to the earthquake epicenter in comparison to ANA, TYH shows anomalous dispersion around two weeks before the earthquake with very weak trend depletion. As mentioned in Section 2, a possible overlapping effect may arise during this outcome due to the presence of another earthquake on November 11, 2016. In a similar study, Asano [39] found significant preseismic impacts for the same earthquake using the same VLF network configuration but a different signal frequency (LF 40 kHz). In that case, the EPZ covers the initial parts of all the propagation paths. Also, as the signal frequency is higher than the frequency of JJI frequency, it can penetrate deeper into the ionosphere and capture more perturbation information. This can allow a possibility of gathering information from different reflection heights, and thus, the change in nighttime amplitude profile behaves in different ways for LF and VLF. This is strong evidence of the seismogenic impression on altitudinal dependency.The overall impression of the VLF anomalies is interesting in nature. It is evident from Figure1 that the earthquake epicenter, receiving locations, and the perturbing zones (The EPZ, CZ, and FZ) interact with each other in a mixed manner. For the EPZ, all the receivers are within the EPZ but the transmitter side of the propagation paths is out of it. In some cases, like ANA, even up to 50% of the paths are outside of the EPZ. Similarly, for TYH, a reasonable portion of such a path has the same features as ANA. Out of all the paths’ 5th FZ, only NSB crosses 20% to 25% of the CZ and KTU marginally touches the CZ. IMZ is a bit far from the CZ and thus as mentioned above, KTU and IMZ show a similar type of moderate change one in the dispersion and another in the trend, respectively. Figure 13 shows the spatial distribution of the signal amplitude as computed by the wave-hop method. It is clear that almost 45% of the total path for ANA is within the skip zone (268 km). This enables us to have minimum information of the skywave and as it is far from the epicenter, ANA does not show any significant changes. On the other hand, having the maximum path length, NSB shows a significant effect. This is due to the fact that NSB is going for multihop reflection and is capable to gather the perturbation information within the CZ by having a reflection point over its path. A detailed statistical analysis of earthquakes in India and its subcontinent region, [30] shows a similar method of using the reflection point over the earthquake perturbed region to get the maximum seismogenic anomalies in the VLF terminator times. The signal strength and path length for TYH and IMZ are close to each other and have similar moderate effects. The situation for KTU is rather different as it stays at one of the minima of the signal amplitude profile. This may be a reason for a weak signature in the dispersion profile.Figure 13
Spatial amplitude profile of JJI signal as a function of distance (in km) as computed from the wave-hop method. The positions of the five receiving locations are marked with vertical lines.The response of the acoustic channel is consistent in all the employed methods. Significant preseismic AGW activity is obtained from VLF fluctuation on November 07-09, and 18-19 for ANA, November 9, 11, 13, and 15-16 for IMZ, November 12 for KTU, November 11, and 21 for NSB, and November 7, 14, and 17-18 for TYH station. In the direct measurement of AGW from the SABER/TIMED temperature profile, the enhancement in Ep associated with the AGW is found to be maximum over the epicenter from November 16 to 19 at 46 km altitude. Also, maximum intensification in wave-like structure (AGW) is obtained from the small-scale fluctuations as observed from GNSS-TEC on November 9 and 18, 2016. The AGW activity on November 07-09 may be connected to the preseismic effect of the second earthquake that took place on November 11. In the direct measurement of higher ionospheric height (F-layer), a noticeable enhancement in TEC is observed on November 12 and 16, 2016. The overall results are consistent in nature with the recent multiparametric findings by [25]. This also corroborates the findings by [47] where the sensitivity of the different parameters of the two major channels in LAIC has been discussed widely. It is well known that AGW can also be generated due to meteorological phenomena, and we investigate this during our study. We find the period November 15 to 20, 2016 is meteorologically quiet (http://www.jma.go.jp/jma/indexe.html). Therefore, the AGW activity during this period is purely seismogenic in nature. The ionospheric TEC highly depends on the solar activity parameters Dst and Kp index. Apart from a moderate storm on November 10, the overall condition was geomagnetically quiet and there is no such contamination in the TEC and VLF profiles.We must mention that the analysis of seismogenic perturbation may give different outcomes for statistical and case study results. In the previous findings (Section1), the employment of different earthquake nucleation zones (e.g., EPZ and CZ) and their overlap with propagation paths’ area (e.g., FZ) are assumed to lead to prominent seismogenic effects. The impression of such seismogenic effect may have different nature in context to the different preparation zone. For example, for VLF radio sounding, the reception center within the EPZ may not show a significant preseismic effect in the signal even if it is within the EPZ. Ghosh et al. [79] show that in the same EPZ, two VLF receiving locations show a sharp difference in the seismogenic anomalies. In statistical and case wise studies, it is obvious that a preparation zone will also tend to diminish the influence of an earthquake as one goes far from the epicenter [80–82]. A statistical analysis inside an EPZ may give a combined outcome of all the earthquakes that their EPZ overlaps with the EPZ of interest, and their occurrence is relatively close in time, having individual weightage depending on the distance of the overlap area from their epicenter, and how close in time they occur. For two simultaneous earthquakes having the same epicenter, the effect of the weakest earthquake may be surpassed by the effect of the strongest due to the concentric preparation zones. Sasmal and Chakrabarti [30] show such an effect by using the concept of the effective magnitude of earthquake when there is overlapping of such preparation zones. For a single earthquake, the scenario may be different for different parameters. As the EPZ, FZ, and CZ are different for the VLF reception stations, significant path and preparation zone dependency is expected to be present in the seismogenic anomalies. This manuscript gives such a mixed example of utilizing different zones where the same VLF signal shows different seismogenic impressions. It is also noticeable that the wavelike structures in night-time signal amplitude give a better result than the conventional night-time fluctuation. The trend and dispersion give a whole average impression of the fluctuation profile in the night-time signal. On the contrary, the wavelet and FFT outcomes gives small-scale fluctuation in the same signal due to overlapping with wavelike structures having the periods of the AGW. Thus, even though VLF signals are being used for the study of two different channels of LAIC, the outcomes may be different in nature.


## 5. Conclusion

This manuscript deals with a multiparametric approach to investigate the acoustic and electromagnetic channels of the LAIC mechanism during the Fukushima earthquake on November 21, 2016. A VLF network in Japan, a GNSS-IGS station (MIZU), and a satellite observation (SABER/TIMED) have been incorporated to check the possible preseismic irregularities in the ionospheric D and F layers and in the stratosphere caused by earthquakes in this study. The VLF signal from the JJI transmitter is recorded at five receiving locations (ANA, IMZ, NSB, TYH, and KTU) to cover a greater geographical area of interest. In the electromagnetic channel, the nighttime signal amplitude shows minimal unusual changes in trend and dispersion. In the same electromagnetic channel, a significant increase in VTEC is observed around one week before the earthquake. These results raise the necessity of using a multiparameter approach for a single-channel observation. The wavelike structures having the periodicity of AGW are found in the VLF nighttime fluctuations using both the FFT and wavelet analysis; whereas, the same VLF data are unable to detect much electromagnetic signature in the lowermost region of the nighttime ionosphere. We also see the altitudinal anisotropy for the first time. Direct evidence of similar waves in the stratospheric region is observed by the SABER/TIMED satellite, and enhancement of the potential energy of such waves is found over the earthquake epicenter. By using proper fitting methods, small-scale fluctuations in VTEC profiles are also found before the earthquake. The overall impression of the study is that for this earthquake, the acoustic channel gives better outcomes in comparison to the electromagnetic channel. However, we found some other studies, other significant parameters such as ultra-low frequency emission (ULF), energetic particle precipitation in the inner radiation belt, enhancement in electron temperature, and electron density using the SWARM satellite showed significant anomalies before the same earthquake. In [25], the anisotropy of different parameters of LAIC has been elaborately mentioned. This manuscript also validates the same as different parameters show different time windows for the preseismic indications. As the VLF outcomes show mixed effects, for a better understanding, VLF signals for similar earthquake locations need to be analyzed and checked with other parameters. This will also enable us to get a better idea of the propagation path and frequency dependency of seismogenic impressions.



---

*Source: 7159204-2023-02-24.xml*