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

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


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][2][3][4][5][6][7][8][9][10][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 preand coseismic anomalies [12][13][14][15][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][21][22], anomalies in the ionospheric total electron content (TEC) [23][24][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][27][28].The statistical studies had been applied by varying the depth and magnitudes of earthquakes [29][30][31][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][52][53][54][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 2 Journal of Sensors 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 magnitude M = 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.

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 Table 1.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 smallscale 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/).
The computation and calibration of TEC are done by the well-known method prescribed by Gopi Seemala by using the biasing error [62][63][64][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 K p, Dst, and a p 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 3 Journal of Sensors 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 fol-lowing equations Here, N b and N e are the start and end times of the nighttime data in seconds.A avg ðtÞ is the average amplitude at a time (t) of 31 nights, and AðtÞ 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   Journal of Sensors 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 oneminute 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.[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][73][74][75][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.gatsinc.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 (T 0 ) 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 (T p ).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,

AGW Computation
where z is the altitude, c p is the specific heat at constant pressure, and g is the acceleration due to gravity.The potential energy (E p ) associated with the AGW for individual temperature profile can be obtained by putting the values of perturbation and background temperature in 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).

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 median X 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 dV TEC, the VTEC is the observed TEC profile, and the VTE C f is the fitted profile, then it can be written as, 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) ..., B n−1 , B n ), 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 (z k ) s as, After getting a small-scale fluctuations or dVTEC, the spectral analysis of dVTEC is 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].

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 Figures 3 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.
By using the FFT and wavelet analysis, we demonstrate the presence of wavelike structures in the nighttime data in Figures 5 and 6. Figure 5  The outcomes of FFT analysis give a mixed impression.For ANA, (Figure 5(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 6 Journal of Sensors The wavelet power spectrum of the ANA station (Figure 6(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).

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.Figure 7 shows that a maximum value of Ep on this day is 11.47 J/kg at an altitude of 44 km.
The temporal variation or time-altitude variation of Ep with the associated AGW activity is also obtained from the nine-dimensional matrix.From Figure 8, 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.
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 timealtitude variation.From Figure 9, it is observed that for   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.

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 Figure 10.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 (Dst min = −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.
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 abovementioned range.The wavelet spectrum as computed from the small-scale fluctuation satisfying the said condition is shown in Figure 11.
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 Figure 12.A comparison of the maximum changes in the context of the temporal weightage of all the methods is tabulated in Table 2.

Discussion
The simultaneous study of ground and space-based methods using two major channels of the LAIC mechanism gives   11 Journal of Sensors 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 (Table 1) show a mixed nature in anomalies.IMZ shows a strong depletion in trend and a medium enhance-ment 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 Figure 1 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 14 Journal of Sensors 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 15 Journal of Sensors 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.
The response of the acoustic channel is consistent in all the employed methods.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 K p 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 (Section 1), 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][81][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  17 Journal of Sensors 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-

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  19 Journal of Sensors 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.

Figure 1 :
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 :
Figure 2: Dst (top panel), K p (middle panel), and a p (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).
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 4 :
Figure 4: Same as Figure 3, but for normalized dispersion.

Figure 5 :Figure 6 Figure 6 :
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.

Figure 6 :Figure 7 :
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.18-11-2016

Figure 8 :Figure 9 :
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.

Figure 10 :
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.

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).

Figure 12 :
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).

Figure 13 :
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.

Table 1 :
Locations of the VLF transmitter and receiver along with the distance from the earthquake epicenter.
Transmitter Receiver Location Distance from the transmitter (km) Distance from the earthquake epicenter (km) JJI (22.1 kHz) 32.04 °N, 130.81 °E ANA 33.89 °N, 134.66 °E [70]g 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.
Significant preseismic AGW activity

Table 2 :
Comparative study of preearthquake effect.