A Multi-Parametric Investigation of Seismogenic Signals: A Case Study of the 2018 Aomori Earthquake (Mw 6.2)

Kousik Nanda^1,2,3, Bhuvnesh Brawar⁴, Masashi Hayakawa^5,6*, Yasuhide Hobara^7,8, Stelios M. Potirakis^9,10, Ajeet K. Maurya¹¹, Uma Pandey¹², Abhirup Datta⁴, Sudipta Sasmal^1,5,6,9*

¹Institute of Astronomy Space and Earth Science, Kolkata, India.
²Department of Space Science and Technology, Adamas University, Kolkata, India.
³Department of Computer Science and Application, Prabhat Kumar College, Contai, India.
⁴Indian Institute of Technology Indore, Indore, India.
⁵Hayakawa Institute of Seismo Electromagnetics, Co., Ltd. (Hi-SEM), UEC Alliance Center, Tokyo, Japan.
⁶QuakeInsight Tokyo, UEC Alliance Center, Tokyo, Japan.
⁷Department of Information and Networking, The University of Electro-Communications, Tokyo, Japan.
⁸Center for Space and Radio Engineering, The University of Electro-Communications, Tokyo, Japan.
⁹Department of Electrical and Electronics Engineering, Ancient Olive Grove Campus, University of West Attica, Egaleo, Greece.
¹⁰Institute for Astronomy, Astrophysics, Space Applications and Remote Sensing, National Observatory of Athens, Metaxa and Vasileos Pavlou, Athens, Greece.
¹¹Department of Physics, Babasaheb Bhimrao Ambedkar University, Vidya Vihar, Lucknow, India.
¹²Department of Earth Sciences, Indian Institute of Technology, Kanpur, India.

DOI: 10.4236/ojer.2026.152002                

Abstract

We present a comprehensive multi-parametric analysis of Lithosphere-Atmosphere-Ionosphere Coupling (LAIC) processes associated with the 24 January 2018 Aomori earthquake ($M~6.2$ ), located off the east coast of northern Honshu, Japan. Thermal observations reveal coherent pre-seismic irregularities in near-surface air temperature, relative humidity, and atmospheric chemical potential (ACP), with notable variations emerging 1 - 2 days prior to the event, followed by co-seismic intensification and rapid post-seismic relaxation. Acoustic channel analysis indicates enhanced atmospheric gravity wave (AGW) activity derived from ERA5 reanalysis, with significant amplification beginning 3 - 5 days before the earthquake, peaking around the event time, and gradually diminishing afterward. Electromagnetic and ionospheric signatures are evident in both lower and upper atmospheric regions. Very-Low-Frequency (VLF/LF) sub-ionospheric observations exhibit statistically significant pre-seismic anomalies, particularly one day prior to the earthquake, while GNSS-based Vertical Total Electron Content (VTEC) shows clear enhancements during the preparatory phase and maximum deviation on the earthquake day. Global Ionospheric Map (GIM) analysis further confirms localized ionospheric perturbations centered around the epicentral region, with a clear spatial decay away from the source. All observations were obtained under geomagnetically quiet conditions, ensuring minimal external influence. Together, these thermal, acoustic, and electromagnetic signatures indicate a consistent sequence of pre-seismic build-up, co-seismic response, and post-seismic recovery, providing a robust multi-channel imprint of the Aomori earthquake and highlighting the importance of integrated multi-parametric approaches for understanding earthquake preparatory processes.

Keywords

Lithosphere-Atmosphere-Ionosphere Coupling, Multi-Parameter Analysis, Ground and Space-Based Study, Pre-, Co-Post-Seismic Impressions

1. Introduction

Earthquakes are complex geophysical phenomena that extend beyond the simple release of accumulated strain energy, involving coupled physical, chemical, and electromagnetic processes across pre-, co-, and post-seismic stages. These processes propagate from the lithosphere into the atmosphere and ionosphere, forming a vertically interconnected system linking the Earth’s interior to near-space environments [1]-[6]. Since the mid-twentieth century, earthquake science has evolved toward a coupled geosystem perspective, where lithospheric disturbances produce measurable atmospheric and ionospheric responses, providing new opportunities for identifying earthquake precursors [7].

The Lithosphere-Atmosphere-Ionosphere Coupling (LAIC) mechanism provides a fundamental framework for understanding these interactions [1] [3]. It describes four primary channels—chemical, thermal, acoustic, and electromagnetic—through which seismogenic disturbances propagate upward. The chemical channel involves ionization and gas emissions, the thermal channel reflects perturbations in heat flux and temperature, the acoustic channel is associated with atmospheric gravity waves (AGWs), and the electromagnetic channel encompasses subsurface currents and ionospheric disturbances. These processes are highly nonlinear, anisotropic, and spatially heterogeneous, necessitating integrated multi-parametric approaches for reliable characterization [2].

Thermal anomalies are among the most consistently observed pre-seismic indicators. Perturbations in the Earth’s radiation budget, along with variations in surface temperature, air temperature, surface latent heat flux (SLHF), outgoing longwave radiation (OLR), relative humidity, and atmospheric chemical potential, have been widely reported prior to major earthquakes [8]-[12]. Satellite observations reveal spatially coherent thermal signatures that often precede seismic events [13]-[16]. These anomalies can induce buoyancy-driven convection, generating AGWs that transport energy to higher atmospheric layers. Such disturbances, produced by stress-induced thermal and pressure variations, are detectable using radar, satellite, GPS, and radio techniques [17]-[22].

The electromagnetic channel spans from subsurface regions to the magnetosphere and exhibits highly complex behavior. Ultra Low Frequency (ULF) emissions near earthquake epicenters are widely recognized as key indicators of pre-seismic activity [23]-[26]. These emissions are often linked to energetic particle bursts in the radiation belts, indicating lithosphere-magnetosphere coupling [27]-[29]. In the lower ionosphere, seismo-electromagnetic disturbances are effectively monitored using sub-ionospheric Very Low Frequency (VLF) radio techniques.

Extensive VLF/LF studies have established strong statistical and physical correlations between ionospheric perturbations and seismic activity. Early investigations demonstrated systematic anomalies in radio signal propagation prior to earthquakes [5] [30] [31], while subsequent works confirmed these findings across different regions and propagation paths [32]-[38]. Recent analyses emphasize quantitative correlations and regional variability, further supporting LAIC processes [39]-[42].

Global Navigation Satellite System (GNSS)-based observations provide additional evidence through variations in Vertical Total Electron Content (VTEC), which frequently exhibit anomalies prior to strong earthquakes [43]-[46]. Such precursory signals have been reported across multiple seismic regions, including Wenchuan, Chile, Nepal, and Tamenglong [47]-[50]. Additional parameters, such as F-layer critical frequency variations, Traveling Ionospheric Disturbances (TIDs), and satellite-based measurements of magnetic field and electron density (e.g., Swarm), further support seismo-ionospheric coupling [2] [51]-[53]. Ground-based observations, including GPS-derived deformation and aerosol anomalies, provide complementary evidence [54] [55].

Despite these advances, pre-seismic anomalies remain highly variable in both space and time, and no single parameter can fully characterize earthquake preparation processes. Consequently, recent research emphasizes multi-parametric approaches that integrate diverse datasets to identify robust precursors. Such frameworks have been successfully applied to major earthquakes worldwide, including the 2021 Samos earthquake (M = 6.9) and events across Japan, Italy, Nepal, and China [2] [40] [56]-[62].

This study introduces an integrated multi-parametric framework that simultaneously analyzes thermal, acoustic, electromagnetic, and ionospheric parameters to identify coherent pre-seismic signatures. By combining ground-based and satellite observations, the approach enables improved characterization of the spatiotemporal evolution of anomalies while explicitly accounting for the nonlinear and event-specific nature of seismogenic processes.

The framework is applied to the Aomori earthquake of 24 January 2018 (19:51 JST; M = 6.2; depth 30 - 40 km), located off the east coast of Aomori Prefecture near Misawa, Japan. This event occurred within the Japan Trench subduction system, where the Pacific Plate subducts beneath the Okhotsk/North American Plate. Its tectonic setting provides an ideal case for investigating lithosphere-atmosphere-ionosphere coupling, with emphasis on identifying pre-seismic anomalies, their temporal evolution, and inter-parameter relationships within a coupled geospace system.

2. Data and Analysis

The earthquake analyzed in this study occurred on 24 January 2018 at approximately 19:51 JST, with a moment magnitude of $M~6.2$ and a focal depth of about 30 - 40 km (https://www.jma.go.jp/jma/index.html). The epicenter was located off the east coast of Aomori Prefecture, near Misawa in northern Honshu, Japan, at geographic coordinates (41.10˚N, 142.43˚E). This event occurred within the tectonically active Japan Trench subduction system, where the Pacific Plate is subducting beneath the Okhotsk (or North American) Plate. The region is characterized by complex plate interactions, high seismicity, and significant stress accumulation, frequently generating moderate- to large-magnitude earthquakes. Owing to its offshore location and relatively shallow focal depth, this event is particularly suitable for investigating lithosphere-atmosphere-ionosphere coupling processes, as stress-induced perturbations can efficiently propagate from the crust to the upper atmospheric and ionospheric layers. The well-defined tectonic framework, combined with the availability of multi-instrument observational data, makes this earthquake an ideal case for examining pre-seismic anomalies and their spatiotemporal evolution within a coupled geophysical system.

(a)

(b)

Figure 1. (a) Location of the earthquake epicenter (red star) and the nearby MIZU IGS station (blue circle) used in this study, along with the estimated EPZ; (b) Great-circle propagation paths of NPM VLF signals to WKN, AMR, and TGN with Fresnel zones, along with the earthquake epicenter and EPZ.

Figure 1(a) shows the location of the earthquake epicenter (red star), along with the nearby MIZU IGS station (blue circle), from which data are used in this study. The blue circle also represents the estimated earthquake preparation zone (EPZ) as computed by the Dobrovolsky formula [63].

VLF/LF observations were conducted by the Hi-SEM network at eleven stations across Japan during the period 2015-2019 [40]. In the present study, data from three representative stations are analyzed: Aomori (AMR), Wakkanai (WKN), and Togane (TGN). Figure 1(b) depicts the great-circle propagation paths of VLF signals transmitted from the NPM 21.4 kHz (Lualualei, Hawaii (21˚25'N, 158˚09'W)) transmitter to three VLF receivers. These are WKN (green diamond, 45.51˚N, 141.93˚E), AMR (cyan diamond, 41.29˚N, 141.21˚E), and TGN (yellow diamond, 35.50˚N, 140.38˚E). The great-circle path lengths are approximately 6007 km for NPM-WKN, 6035 km for NPM-AMR, and 6112 km for NPM-TGN. The corresponding Fresnel zones (similar color-coded curves) represent the regions of maximum wave sensitivity within the Earth-Ionsophere Waveguide EIWG, where perturbations in the lower ionosphere can significantly influence signal amplitude and phase. The earthquake epicenter is marked by the red star. The dashed circle denotes the EPZ. Among the three receivers, AMR is the closest to the earthquake epicenter, WKN lies near the boundary of the EPZ (borderline), while TGN is located outside the EPZ. This spatial configuration allows comparative assessment of ionospheric perturbations as a function of distance from the seismic preparation zone. Importantly, substantial portions of the Fresnel zones intersect the EPZ for all three propagation paths, indicating that the VLF signals are likely modulated by localized ionospheric disturbances associated with the seismic activity.

The solar wind parameters and geomagnetic indices derived from the NASA-OMNIWeB database https://omniweb.gsfc.nasa.gov/ were examined over a 23-day interval spanning 15 days before and 7 days after the earthquake of 24 January 2018 (DOY ~24) (Figure 2). Around the earthquake, geomagnetic conditions remain relatively quiet, with Dst varying between ~−20 to +10 nT, indicating the absence of storm-time activity. The Kp index mostly stays below 3, and the Ap index remains low (<20 nT, with daily Ap ~5 - 10 nT), confirming quiet to mildly disturbed conditions. The F10.7 solar flux is stable around ~66 - 69 sfu, suggesting steady solar EUV input without significant enhancements. The IMF Bz component fluctuates within ~2 - 8 nT without prolonged southward excursions. Although minor increases in Kp (~4) and Ap (~25 - 30 nT) are seen just after the earthquake, no major geomagnetic storm conditions are present, indicating that the observed seismogenic impressions in various channels are unlikely to be driven by external solar-geomagnetic forcing.

Figure 2. Temporal variations of geomagnetic and solar parameters during 09-31 January 2018 indicate predominantly quiet conditions around the earthquake day (24 January), with no significant storm-time disturbances.

2.1. Thermal Parameter Analysis

To examine possible anomalies in surface temperature, relative humidity (RH), and atmospheric chemical potential (ACP), reanalysis data were obtained from the NOAA Physical Sciences Laboratory (https://psl.noaa.gov/data/gridded/data.ncep.reanalysis.html). The analysis was carried out over the region 39˚N - 43˚N and 140˚E - 145˚E, representing an area of approximately 5˚ around the earthquake epicenter. Since the NOAA reanalysis is provided as global fields, the required variables were first downloaded and then subset for the study region. The ACP ($\text{Δ}U$ ) was then estimated from air temperature and RH according to [7]:

$\Delta U=5.8\times {10}^{-10}{\left(20T_g+5463\right)}^2\ln \left(\frac{100}{H}\right),$ (1)

where $T_g$ denotes air temperature and $H$ denotes relative humidity. Spatiotemporal maps were generated using the daily mean values of each parameter over the selected period surrounding the earthquake.

2.2. AGW Analysis

To investigate the potential excitation of AGWs in relation to the seismic event, vertical temperature profiles were retrieved from the ERA5 reanalysis database available via the Copernicus Climate Data Store (https://cds.climate.copernicus.eu). The temperature data span 137 geopotential levels, extending from level 137 (approximately 10 m above the surface) up to level 1 (about 79.30 km altitude), thereby covering the troposphere, stratosphere, and mesosphere. The analysis period includes 15 consecutive days, comprising seven days before the earthquake, the day of occurrence, and seven days afterward, with a temporal resolution of one hour.

AGWs in the middle atmosphere are typically characterized by relatively long vertical wavelengths, often reaching up to ~30 km [64]. However, employing such a broad filtering range would encompass nearly the entire mesosphere (approximately 45 - 80 km) and may introduce contamination from large-scale background structures, such as temperature inversions near the tropopause and stratopause. To mitigate this, a Chebyshev Type-I digital band-stop filter with a stopband corresponding to vertical wavelengths of 2 - 10 km (the typical vertical scales of stratospheric AGWs) and a ripple of 3 dB was applied to the ERA5 temperature profiles $T$ [19] [21]. This filtering step removes gravity-wave-scale components and yields the background temperature field $\overline{T}$ . The perturbation component associated with AGWs was then obtained as:

$T^{\prime }=T-\overline{T}.$ (2)

The atmospheric stability was quantified using the Brunt-Väisälä frequency ($N$ ), which describes the oscillatory response of an air parcel in a stably stratified medium:

$N^2=\frac{g}{\overline{T}}\left(\frac{\text{d}\overline{T}}{\text{d}z}+\frac{g}{c_p}\right),$ (3)

where $z$ denotes altitude (km), $c_p=1.005\text{kJ}\cdot {\text{kg}}^{-1}\cdot {\text{K}}^{-1}$ is the specific heat at constant pressure, and $g=0.0098\text{km}\cdot {\text{s}}^{-2}$ is the gravitational acceleration. The

vertical gradient term $\frac{\text{d}\overline{T}}{\text{d}z}$ was evaluated at a resolution of 100 m using the

background temperature profile.

The gravity wave potential energy per unit mass ($E_p$ ) was subsequently estimated as:

$E_p=\frac{1}{2}{\left(\frac{g}{N}\right)}^2{\left(\frac{T^{\prime }}{\overline{T}}\right)}_{\text{layer}}^2,$ (4)

where the variance term is altitude-dependent and defined as:

${\left(\frac{T^{\prime }}{\overline{T}}\right)}^2=\frac{1}{z_{\text{max}}-z_{\text{min}}}{\displaystyle {\int }_{z_{\text{min}}}^{z_{\text{max}}}{\left(\frac{T^{\prime }}{\overline{T}}\right)}^2\text{d}z}.$ (5)

here, $z_{\text{min}}$ and $z_{\text{max}}$ represent the lower and upper limits of the selected altitude interval. In this study, the variance was computed within discrete 100 m layers to capture fine-scale vertical structures. The temporal evolution of $E_p$ was examined within the altitude range of 35 - 50 km above the earthquake epicenter, corresponding to the upper stratosphere-lower mesosphere region where decreasing atmospheric density enhances gravity wave amplitudes and improves detectability. Furthermore, spatial distributions were constructed by mapping the daily mean $E_p$ values within ±10˚ of the epicenter at an altitude of around 46.3 km. This specific altitude corresponds to the level at which gravity wave potential energy exhibited a pronounced and coherent maximum during the study period, thereby providing a representative level for interpretation, consistent with previous findings [2] [21].

2.3. VLF Analysis

To examine the VLF fluctuations, we use the conventional Nighttime Fluctuation Method (NFM), which is one of the most widely adopted statistical techniques for analyzing VLF signal perturbations [40]. The methodology consists of three principal stages.

First, the raw VLF amplitude data $A\left(t\right)$ (in dB) are used to compute the residual nighttime signal $dA\left(t\right)$ by subtracting the mean nighttime amplitude $\langle A\left(t\right)\rangle$ calculated over a 31-day window (i.e., ±15 days around the day of interest, including the day itself). Thus,

$dA\left(t\right)=A\left(t\right)-\langle A\left(t\right)\rangle .$ (6)

In this study, the nighttime interval is defined as 18:00 - 02:00 UT.

Second, three daily parameters are derived from $dA\left(t\right)$ to characterize VLF signal variability: the trend (TR), dispersion (DP), and nighttime fluctuation (NF). These are defined as:

$\text{TR}=\frac{1}{N_e-N_s}{\displaystyle \sum }_{t=N_s}^{N_e}dA\left(t\right),$ (7)

where TR represents the mean value of $dA\left(t\right)$ over the selected nighttime interval, and $N_s$ and $N_e$ denote the starting and ending time indices, respectively.

The dispersion (DP), corresponding to the standard deviation of $dA\left(t\right)$ , is given by:

$\text{DP}=\sqrt{\frac{1}{N_e-N_s}{\displaystyle \sum }_{t=N_s}^{N_e}{\left(dA\left(t\right)-\text{TR}\right)}^2}.$ (8)

The nighttime fluctuation (NF), representing the signal power, is defined as:

$\text{NF}={\displaystyle \sum }_{t=N_s}^{N_e}{\left(dA\left(t\right)\right)}^2.$ (9)

Third, the daily time series of TR, DP, and NF are normalized to obtain ${\text{TR}}^{\text{*}}$ , ${\text{DP}}^{\text{*}}$ , and ${\text{NF}}^{\text{*}}$ using:

$X^{\text{*}}=\frac{X-{\langle X\rangle }_{\pm 15}}{{\sigma }_{\pm 15}},$ (10)

where ${\langle X\rangle }_{\pm 15}$ and ${\sigma }_{\pm 15}$ represent the mean and standard deviation computed

over the same ±15 day window. Anomalies exceeding the $\pm 2\sigma$ threshold in these normalized series are considered potential earthquake-related disturbances.

It is important to note that TR and DP are generally independent quantities, with TR often regarded as the primary indicator, whereas NF reflects a combined contribution of both trend and dispersion. The above implementation of NFM, based on a symmetric ±15 day window, is standard in retrospective (a posteriori) analyses [2] [38] [40]. However, such an approach inherently incorporates future information. For real-time applications or earthquake forecasting systems, the statistical baseline should instead be computed using only prior data, typically 30 days preceding the day of interest.

2.4. GPS-TEC Analysis

The ionospheric TEC quantifies the integrated electron density along the propagation path between a GNSS and a ground-based receiver. In this study, TEC values were derived from dual-frequency GPS observations using the Gopi Seemala GPS-TEC software (https://seemala.blogspot.com), which has been extensively utilized and validated in ionospheric research [65] [66]. The software computes Slant TEC (STEC) through the geometry-free linear combination of dual-frequency pseudorange and carrier-phase measurements. To enhance accuracy, carrier-phase observations are employed to level the code-derived TEC, thereby reducing measurement noise.

Satellite differential code biases (DCBs) are corrected using International GNSS Service (IGS) DCB products, while receiver DCBs are estimated internally during the processing stage. The resulting STEC values are expressed in TEC units (TECU), where 1 TECU = 10¹⁶ electrons/m². Subsequently, Vertical TEC (VTEC) is obtained by mapping STEC to the vertical direction using a single-layer ionospheric model with a fixed shell height of 350 km, corresponding to the typical altitude of the ionospheric F-region peak [43].

The VTEC values are calculated at the ionospheric pierce points (IPPs) associated with each satellite-receiver line of sight, enabling spatial characterization of ionospheric variability. To ensure data quality and minimize uncertainties arising from multipath effects and low-elevation signal distortions, only observations with satellite elevation angles exceeding 30˚ are considered. These processing procedures and bias corrections ensure the reliability and robustness of the derived TEC measurements for subsequent statistical and seismo-ionospheric analyses.

To identify TEC anomalies, a statistical approach based on the median and interquartile range (IQR) was employed. For each universal time (UT), the median ($X$ ) of VTEC values over a sliding window of 15 days (7 days before, the day of the event, and 7 days after) was computed along with the corresponding IQR. The upper bound (UB) and lower bound (LB) were then defined as:

$\text{UB}=X+\text{IQR},$ (11)

$\text{LB}=X-\text{IQR}.$ (12)

This method is consistent with the approach of [44], where for a normally distributed dataset with mean ($\mu$ ) and standard deviation ($\sigma$ ), the median and IQR correspond approximately to $\mu$ and $1.34\sigma$ , respectively [67]. An anomaly is identified when the VTEC exceeds the upper bound or falls below the lower bound with a confidence level of approximately 80% - 85%. Accordingly, positive and negative anomalies are defined as:

$\text{anomaly}=\text{VTEC}\ge \text{UB},$ (13)

$\text{anomaly}=\text{VTEC}\le \text{LB}.$ (14)

This procedure enables robust detection of statistically significant enhancements and depletions in TEC relative to the background ionospheric state [68].

The GIM analysis in this study is based on the global ionospheric map product available through the JPL GUARDIAN system (https://guardian.jpl.nasa.gov/analysis/globalIonoMovie/index.html). These maps are generated using GPS and Galileo observations within a global differential GNSS framework, with a temporal resolution of 15 minutes and a spatial resolution of 1˚ × 1˚. The analysis was carried out over the region 39 - 43˚N and 140 - 145˚E during 17-31 January 2018. As GIM data represent a model-assimilated global ionospheric state, the interpretation is based on relative spatio-temporal variations within the study period rather than explicit subtraction from a multi-year climatological baseline. Mean VTEC values were obtained by averaging the GIM grid values over the selected study region, enabling characterization of regional ionospheric evolution around the earthquake period.

3. Main Results

3.1. Thermal Irregularities

Figure 3(a) illustrates the daily air temperature anomalies over the region 39 - 43˚N and 140 - 145˚E during 17-31 January, with values ranging approximately from −6 K to +6 K. The anomaly baseline was defined as a multi-year climatological mean computed at each grid point using corresponding daily averages from 2017 and 2019, thereby accounting for both spatial and temporal variability. At the epicentral location (~41.1˚N, 142.4˚E), the anomaly is nearly neutral on 17 January (~0 K), followed by weak cooling on 18 January (~−1 K) and slight warming on 19 January (~+0.5 K). Variations remain small up to 21 January (~+1 K), indicating a weak warming phase, before increasing to ~+2 K on 22 January and then dropping to ~−1 K on 23 January, corresponding to a rapid change of about 3 K within one day. On 24 January (earthquake day), a clear intensification is observed, with the anomaly rising sharply to ~+4 K at the central location, marking the onset of a strong warming event; this peaks on 25 January (~+5 K) and remains elevated on 26 January (~+4 K). The maximum anomaly reaches about +6 K in the southwestern sector (39 - 40˚N, 140 - 141˚E), while the northeastern region shows relatively lower values (~+3 K), producing a pronounced latitudinal gradient of approximately 2 K per degree. Subsequently, after 27 January (~+2 K), a rapid cooling phase occurs, with anomalies decreasing to ~−1 K on 28 January and reaching minimum values near −4 K in the northeast, followed by stabilization during 29-31 January within the range of 0 to +0.5 K with weak spatial variation. Overall, the results indicate a distinct and abrupt warming anomaly initiated on the earthquake day, followed by a rapid post-event relaxation.

(a) (b)

(c)

Figure 3. (a) Daily air temperature anomalies over 39 - 43˚N and 140 - 145˚E during 17-31 January 2018 show a warming phase beginning 2 days prior to the earthquake, peaking on the event day (24 January), and persisting for 2 - 3 days afterward before relaxing; (b) Daily relative humidity (RH) anomalies for the same temporal and spatial span as like Figure 3(a). It clearly shows humid conditions much before the earthquake, a decrease on 23 January, renewed moistening on the earthquake day (24 January), and a subsequent drying post-seismic phase; (c) Daily ACP anomalies for the same temporal and spatial span as in Figure 3(a) show a distinct enhancement immediately prior to the earthquake (23 January), followed by a reduction on the earthquake day (24 January), and a subsequent post-seismic increase, indicating a transient modulation.

Figure 3(b) shows daily relative humidity (RH) anomalies (in %) over the region 39 - 43˚N and 140 - 145˚E during 17-31 January, with values ranging approximately from −20% to +20%. At the epicentral location (~41.1˚N, 142.4˚E), conditions begin with pronounced negative anomalies on 17 January (~−5% to −7%), intensifying further on 18 January (~−8% to −10%). A sharp transition to positive anomalies occurs on 19 January (~+10% to +15%), followed by sustained humid conditions during 20-22 January (~+5% to +12%). On 24 January (earthquake day), RH remains elevated at approximately +5 to +8%, indicating a persistent moist anomaly over the region. Immediately afterward, however, a marked drying phase develops: by 25 January, values decrease to ~0% to −2%, followed by further reduction to ~−6% to −8% on 26 January. The minimum anomalies reach about −15% to −20%, particularly over the southern sector (39 - 40˚N), while relatively higher values (~−3% to 0%) persist in the northern region, producing a latitudinal gradient of approximately 2% - 3% per degree. During 27-28 January, dry conditions dominate (~−8% to −12%), followed by a gradual recovery, with anomalies reaching ~−2% to 0% on 29 January and transitioning back to moist conditions during 30-31 January (~+2 to +8%), with stronger positive anomalies in the northeastern sector (~+10 to +15%). Overall, the observed pattern indicates a humid phase around the earthquake day, followed by a pronounced post-event drying and subsequent recovery.

Figure 3(c) presents daily variations of atmospheric chemical potential (ACP, in eV) over the region 39 - 43˚N and 140 - 145˚E during 17-31 January, with values ranging approximately from −0.02 to +0.02 eV. At the epicentral location (~41.1˚N, 142.4˚E), ACP values are initially positive on 17 January (~+0.004 to +0.006 eV) and show a slight increase on 18 January (~+0.006 to +0.008 eV). A sharp transition to negative values occurs on 19 January (~−0.004 to −0.006 eV), followed by more pronounced negative anomalies on 20 January (~−0.006 to −0.008 eV). On 21 January, the central region remains weakly negative (~−0.003 to −0.005 eV), while 22 January exhibits mixed conditions with a strong zonal gradient, characterized by negative anomalies in the western sector and positive anomalies in the eastern part. On 24 January (earthquake day), ACP values near the epicenter are close to neutral to slightly negative (~−0.002 to 0 eV), indicating a transitional state. Immediately afterward, a clear positive phase develops, with values of ~+0.002 to +0.004 eV on 25 January, further intensifying during 26-27 January (~+0.004 to +0.006 eV). The southern sector (39 - 40˚N) exhibits relatively higher values (~+0.006 to +0.008 eV), producing a latitudinal gradient of approximately 0.001 - 0.002 eV per degree. From 28 January onward, ACP values begin to decrease, showing near-neutral to weakly positive conditions (~0 to +0.002 eV), followed by a pronounced negative phase during 30-31 January (~−0.006 to −0.01 eV), particularly over the eastern region. Overall, the results indicate a transition from positive to negative anomalies prior to the earthquake, near-neutral conditions on the event day, followed by post-event positive enhancement and subsequent decline.

3.2. Acoustic Irregularities

Figure 4 presents the temporal-altitudinal evolution of AGW potential energy per unit mass ($E_p$ ) derived from ERA5 reanalysis over the altitude range ~35 - 50 km during 17-31 January, with the vertical dashed line marking the earthquake occurrence on 24 January 2018. The AGW activity is predominantly confined to the middle atmospheric region (approximately 42 - 48 km), where the majority of energy is concentrated. During the pre-seismic phase (17-22 January), moderate wave activity is observed, with $E_p$ values generally in the range of ~10 - 25 J/kg and localized enhancements up to ~28 J/kg around 22-23 January near 44 - 46 km altitude. A clear intensification begins just prior to the earthquake, with $E_p$ increasing to ~26 - 30 J/kg on 23-24 January. On the earthquake day, enhanced AGW activity is evident, particularly within the 43 - 47 km altitude band. The most significant feature occurs immediately after the event, where a strong post-seismic amplification is observed on 25 January, with peak $E_p$ values exceeding ~35 - 38 J/kg, indicating substantial energy enhancement. This increase exhibits a well-defined vertical structure, with a dominant maximum around 43 - 45 km and secondary extensions up to ~48 km, suggesting efficient upward propagation of wave energy. Following this peak, the AGW activity gradually diminishes during 26-27 January (~15 - 25 J/kg) and further decreases to background levels (~2 - 10 J/kg) after 28 January. Below ~40 km altitude, $E_p$ values remain consistently low (<10 J/kg), indicating that the dominant wave activity is concentrated in the upper stratosphere–lower mesosphere region. Overall, the results demonstrate a clear enhancement of AGW energy beginning shortly before the earthquake, peaking immediately after the event, and subsequently decaying, consistent with the generation and vertical propagation of seismically induced atmospheric disturbances within the lithosphere-atmosphere-ionosphere coupling framework.

Figure 4. Temporal-altitudinal evolution of AGW potential energy ($E_p$ ) derived from ERA5 over 35 - 50 km during 17-31 January 2018, with the earthquake day (24 January) marked by a dashed line. Enhanced AGW activity is observed in the 42 - 48 km region, showing pre-seismic intensification and a pronounced post-seismic peak (>35 J/kg) followed by gradual decay.

3.3. Electro-Magnetic Irregularities

3.3.1. VLF Signal Irregularities

Figures 5(a)-(c) illustrate the results of the Nighttime Fluctuation Method (NFM) where the top panel shows the normalized trend (${\text{TR}}^{\text{*}}$ ), representing the nighttime average amplitude; the middle panel shows the normalized dispersion (${\text{DP}}^{\text{*}}$ ); and the bottom panel presents the normalized nighttime fluctuation (${\text{NF}}^{\text{*}}$ ).

Focusing on the AMR station (Figure 5(a)), which is closest to the earthquake epicenter, two pronounced anomalies are observed in the most significant parameter, ${\text{TR}}^{\text{*}}$ , with values dropping below the $-2\sigma$ threshold approximately 5 days and 1 day prior to the earthquake. A comparison with the geomagnetic Dst index (Figure 2) indicates that the earlier depletion (around 5 days before the event) is associated with a weak geomagnetic storm occurring on 19 January, while the subsequent anomaly around 4 days prior can be interpreted as a residual or recovery-phase effect, consistent with the known persistence of geomagnetic disturbances. In contrast, the depletion observed one day before the earthquake is not correlated with geomagnetic activity and is therefore considered to be of seismogenic origin. Although the magnitude of this anomaly is moderate, likely due to the relatively small earthquake magnitude ($M=6.2$ ), it is accompanied by a simultaneous enhancement in ${\text{DP}}^{\text{*}}$ , which is a characteristic signature of VLF/LF precursory behavior.

(a)

(b)

(c)

Figure 5. (a) Time series of normalized ${\text{TR}}^{\text{*}}$ , ${\text{DP}}^{\text{*}}$ and ${\text{NF}}^{\text{*}}$ at the AMR station showing anomalies below the $-2\sigma$ threshold, with the pre-seismic anomaly (1 day before the event) identified as likely of seismogenic origin; (b) Same as Figure 5(a) for WKN station; (c) Same as Figure 5(a) for TGN station.

For the WKN station (Figure 5(b)), located at the highest latitude among the three stations, the ${\text{TR}}^{\text{*}}$ parameter exhibits a weaker depletion one day prior to the earthquake, which can still be interpreted as a seismogenic effect, albeit less pronounced than at AMR. However, more significant depletions observed 5 - 6 days before the earthquake are clearly associated with the geomagnetic disturbance. It is well established that high-latitude stations, such as WKN, are more sensitive to geomagnetic activity, and therefore storm-related signatures are more prominently detected at this location.

In contrast, the TGN station (Figure 5(c)), situated at a lower latitude and at a greater distance from the epicenter, shows a distinct depletion in ${\text{TR}}^{\text{*}}$ one day before the earthquake, which is attributed to seismogenic ionospheric perturbations. Notably, geomagnetic storm effects are nearly absent at this station, further supporting the non-geomagnetic origin of the observed anomaly. Despite its greater epicentral distance compared to WKN, the presence of a clear anomaly at TGN suggests that the spatial extent of the seismogenic ionospheric disturbance extends significantly toward lower latitudes. This behavior is consistent with previous observations of major seismic events, such as the 2011 Tohoku earthquake [37] [38], in which ionospheric perturbations were predominantly observed over lower-latitude regions.

3.3.2. GPS-TEC Irregularities

Figure 6(a) illustrates the temporal variation of VTEC and its anomaly over the period 17-31 January, with the vertical dashed line marking the EQ occurrence on 24 January 2018. Panel (a) shows the observed VTEC (black curve) along with its corresponding mean (blue curve), exhibiting a clear diurnal pattern with amplitudes typically ranging between ~5 and 12 TECU. This regular oscillation reflects normal ionospheric behavior driven by solar forcing; however, noticeable deviations from the mean are observed during specific intervals. In particular, around 23-24 January, VTEC shows enhanced fluctuations, with peak values reaching ~13 - 14 TECU, exceeding the usual diurnal range. Panel (b) presents the corresponding VTEC anomaly, highlighting departures from the mean behavior. During the pre-seismic period, moderate anomalies are observed, generally within ~0 - 2 TECU, with occasional enhancements up to ~2.5 TECU around 19-21 January. A significant increase is observed on 23 January, where the anomaly reaches ~3.5 - 4 TECU, followed by a pronounced peak on the EQ day (24 January), exceeding ~4.5 TECU, representing the maximum deviation during the entire interval. After the EQ, the anomaly rapidly decreases, returning to lower levels (<1 TECU) during 25-26 January, with only sporadic fluctuations thereafter. Overall, the results indicate a clear intensification of ionospheric variability immediately prior to and during the EQ, followed by a rapid relaxation, suggesting a possible seismo-ionospheric disturbance superimposed on the background diurnal variation.

(a)

(b)

(c)

Figure 6. (a) Temporal variation of the Vertical Total Electron Content (VTEC) recorded at the MIZU00JPN station from Day of Year (DoY) 17 to 31. The (top) panel shows the observed VTEC profile (black curve) along with the upper (UB, red curve) and lower (LB, green curve) bounds. The earthquake day (DoY 24) is marked with a yellow line. The (bottom) panel presents the corresponding TEC fluctuations, expressed as deviations from the background values; (b) Daily GIM VTEC over 39 - 43˚N and 140 - 145˚E shows localized enhancement near the epicenter before and after the earthquake, with a reduction on the earthquake day and subsequent recovery to background levels; (c) Mean VTEC variation with distance from the earthquake epicenter on 23-25 January 2018 shows a decreasing trend with distance, with pre-seismic enhancement, co-seismic reduction, and post-seismic increase.

Figure 6(b) shows daily GIM VTEC over 39 - 43˚N and 140 - 145˚E during 17-31 January 2018, with the black dot marking the EQ epicenter. A clear latitudinal gradient is observed, with higher VTEC (~14 - 20 TECU) in the southern sector (39 - 41˚N) and lower values (~2 - 8 TECU) toward the north (>42˚N). During 17-22 January, VTEC remains relatively stable with moderate variability. A noticeable enhancement appears on 23 January near the epicenter (~15 - 18 TECU), followed by sudden decrease on the EQ day (24 January). A stronger enhancement is observed during 25-26 January, with peak values exceeding ~18 - 20 TECU, mainly in the southwestern region, along with a sharper latitudinal gradient. After 27 January, VTEC gradually decreases, returning to background levels (~10 - 15 TECU) by 29-31 January. Overall, the results indicate localized VTEC enhancement around the EQ period, suggesting possible seismo-ionospheric perturbations superimposed on the background structure.

The variation of mean VTEC (TECU) with distance from the EQ epicenter for 23 (blue), 24 (orange), and 25 (green) January 2018 is shown in Figure 6(c). The mean VTEC values presented in Figure 6(c) were obtained by calculating the daily mean GIM value around the earthquake epicenter. The study area was divided into 10 km radial bins, and all VTEC observations within each bin were averaged to produce a representative mean value for that distance. This procedure was applied separately for 23, 24, and 25 January 2018, enabling the assessment of pre-, co-, and post-seismic ionospheric variations. A clear decreasing trend of VTEC with increasing distance is observed for all three days, indicating that the ionospheric response is strongest near the epicentral region. On 23 January (pre-EQ), VTEC values are relatively elevated near the epicenter (~14.5 - 14.8 TECU within ~10 - 30 km) and gradually decrease to ~12.7 TECU at ~95 km. On the EQ day (24 January), an overall reduction in VTEC is evident, with values decreasing from ~13.8 - 13.9 TECU near the epicenter to ~12.4 TECU at larger distances, suggesting a temporary suppression of ionospheric electron density. In contrast, on 25 January (post-EQ), a clear enhancement is observed, with peak VTEC reaching ~15.6 TECU near ~15 - 25 km, followed by a gradual decline to ~13.5 TECU at ~95 km. This behavior is consistent with the spatial pattern observed in GIM maps, where enhanced TEC is concentrated near the epicenter and diminishes outward. The systematic radial decrease, together with pre-EQ enhancement, EQ-day depletion, and post-EQ amplification, suggests a localized ionospheric perturbation centered around the epicentral region, likely associated with LAIC-related processes.

4. Discussions

The results presented in this study provide a coherent and quantitatively consistent multi-parametric perspective of lithosphere-atmosphere-ionosphere coupling processes associated with the 24 January 2018 Aomori earthquake. The thermal parameters (temperature, RH, and ACP) exhibit systematic and measurable variations around the event. Surface air temperature shows a sharp increase of approximately 4 - 5 K on the earthquake day relative to the background level (~0 - 1 K), corresponding to an enhancement of nearly 300% - 400%. This is followed by a rapid post-seismic decrease to ~−1 to −4 K, indicating a strong relaxation phase. In contrast, RH displays an inverse response, with a pre-seismic transition from negative anomalies (~−10%) to positive values (~+10% - 15%), followed by a drop of nearly 15% - 20% after the earthquake. ACP variations further support this coupling, showing a pre-seismic decrease from ~+0.006 eV to ~−0.006 eV (a relative change exceeding 150%), near-neutral conditions on the earthquake day, and a post-seismic increase of ~0.004 - 0.006 eV before declining again. These combined thermal signatures indicate significant modulation of atmospheric thermodynamic conditions associated with seismic processes.

The acoustic channel, represented by AGW activity, exhibits a clear quantitative enhancement in $E_p$ prior to and following the earthquake. During the pre-seismic phase, $E_p$ values increase from background levels of ~10 - 15 J/kg to ~26 - 30 J/kg (an increase of ~100% - 150%), followed by a peak exceeding ~35 - 38 J/kg immediately after the earthquake, corresponding to an overall enhancement of nearly 200% - 250% relative to background conditions. The vertical confinement of this enhancement within 42 - 48 km altitude indicates efficient upward propagation of energy. The subsequent decrease to ~2 - 10 J/kg after 28 January reflects a return to background levels, confirming the transient nature of the disturbance.

Electromagnetic and ionospheric observations further demonstrate strong quantitative signatures. VLF/LF measurements show pre-seismic ${\text{TR}}^{\text{*}}$ depletions exceeding the $-2\sigma$ threshold, indicating statistically significant anomalies. GNSS-based VTEC analysis reveals anomaly amplitudes increasing from background levels of ~1 - 2 TECU to ~4.5 TECU on the earthquake day, representing an enhancement of over 100% - 150%. Similarly, GIM observations indicate localized TEC increases from ~14 - 15 TECU to ~18 - 20 TECU (approximately 20% - 30% enhancement) in the epicentral region. Radial analysis further confirms a spatial decay of ~2 - 3 TECU (about 15% - 20%) within ~100 km from the epicenter, demonstrating that the perturbation is strongly localized.

A comparative assessment of different parameters reveals a consistent and quantifiable sequence of pre-, co-, and post-seismic responses. Pre-seismic signatures are characterized by moderate thermal variations (~1 - 2 K), RH increase of ~10% - 15%, ACP reversal (~150% change), AGW enhancement of ~100%, and ionospheric anomalies exceeding ~2 - 3 TECU. Co-seismic effects are marked by peak thermal anomalies (~4 - 5 K), maximum VTEC deviations (~4.5 TECU), and transitional ACP behavior. Post-seismic responses include strong AGW amplification (up to ~250%), RH depletion of ~15% - 20%, and gradual relaxation of ionospheric parameters to background levels. The temporal alignment and quantitative consistency among these parameters strongly support a coupled geophysical response rather than isolated variations.

A consistent significance assessment across all parameters confirms the robustness of the observed anomalies. Thermal parameters show clear exceedance of their respective thresholds, with temperature anomalies of ~4 - 5 K surpassing $\pm 2\sigma$ (~±2.4 K), RH variations of ~15% - 20% exceeding $\pm 2\sigma$ (~±8%), and ACP changes beyond $\pm 2\sigma$ (~±0.0036 eV). Similarly, the acoustic channel exhibits strong enhancement, with AGW energy increasing to ~35 - 38 J/kg, well above the $\pm 2\sigma$ threshold (~±9 J/kg). The ionospheric response is also significant, as VTEC anomalies reach ~4.5 TECU, exceeding $\pm 2\sigma$ (~±3 TECU), while GIM-based VTEC enhancements up to ~18 - 20 TECU surpass the regional $\pm 2\sigma$ threshold (~±6 TECU). Collectively, these statistically significant perturbations across thermal, acoustic, and ionospheric channels demonstrate coherent multi-parameter variability associated with the earthquake, supporting a coupled lithosphere-atmosphere-ionosphere response.

Overall, the multi-parametric analysis demonstrates that integrating thermal, acoustic, and electromagnetic observations provides a more robust and physically consistent framework for understanding earthquake-related processes. The observed magnitudes of variation across domains highlight the efficiency of energy transfer within the LAIC system and underscore the importance of combined-parameter analysis for identifying the complex physics underlying earthquake preparation processes.

In the ~15 days preceding the 24 January 2018 Aomori-region earthquake, environmental conditions were governed by typical winter variability, with no typhoon or tropical cyclone activity affecting northern Japan. The only notable meteorological event was a short-lived winter snowstorm on 22-23 January, associated with a passing extratropical low-pressure system, moderate pressure fluctuations (~5 - 10 hPa), and cold-air advection—conditions that are climatologically common for Aomori in January (https://www.jma.go.jp/jma/indexe.html). Additionally, a minor geomagnetic disturbance around 19 January (likely substorm-level activity) remained within low-to-moderate geomagnetic indices (e.g., $K_p\le 4$ ), insufficient to produce strong ionospheric or lithospheric coupling effects at regional scales (https://wdc.kugi.kyoto-u.ac.jp/). Taken together, both meteorological and geomagnetic conditions during this period fall within background seasonal norms, lacking the magnitude, persistence, or spatial coherence required to account for sustained anomalies, thereby supporting a seismogenic origin for the observed signals.

This study is based on a retrospective analysis of a single earthquake event. However, the reported anomalies are not evaluated in isolation; they are derived relative to corresponding non-earthquake periods under geomagnetically and meteorologically quiet conditions, thereby incorporating a robust reference baseline. Furthermore, the observed variations across all parameters are consistent with their known long-term statistical behavior under non-seismic conditions, as established in previous studies. The multi-parameter anomalies also exhibit strong temporal alignment and internal consistency across independent datasets, supporting their close association with earthquake-related processes. Although the precise causal mechanism governing these lithosphere-atmosphere-ionosphere interactions remains to be fully established, the coherence and quantitative agreement among the observed signatures provide robust observational evidence of a coupled geophysical response.

5. Conclusions

The results obtained in this study demonstrate clear and consistent multi-parametric signatures associated with the seismic event. A notable and potentially novel feature of the present analysis is the asymmetric behavior of the observed anomalies, characterized by pronounced enhancements during the pre- and post-seismic phases, while the co-seismic period exhibits a relative depletion across multiple parameters. This behavior contrasts with earlier studies (see Introduction), where co-seismic enhancements or mixed responses were more commonly reported. In the present case, the suppression of anomalies on the earthquake day, despite the absence of significant geomagnetic or meteorological disturbances, suggests a transient reorganization of the LAIC system.

Interestingly, the presence of both pre- and post-seismic anomalies has been reported in several previous studies. Statistical and case-based analyses indicate that ionospheric perturbations can occur several days before and also persist for a few days after an earthquake, typically within a window of ~5 days prior and up to ~3 - 5 days after the event [69] [70]. For example, multi-event studies have identified TEC anomalies appearing 4 - 10 days before and 3 - 5 days after earthquakes, supporting the existence of both preparatory and relaxation phases of the ionospheric response [71] [72]. Thus, post-seismic ionospheric disturbances can also be interpreted as either residual effects or secondary adjustments of the coupled system following the main shock. However, the nearly symmetric and comparable magnitude of pre- and post-seismic responses observed in the present study is relatively uncommon. This suggests that the processes governing stress accumulation and post-rupture relaxation may follow similar physical pathways, possibly involving comparable energy transfer mechanisms through the LAIC system. Such symmetry implies that post-seismic anomalies are not merely residual signatures but may represent an integral part of the overall coupling dynamics. These findings highlight the importance of considering both pre- and post-seismic intervals within a unified framework and suggest a more complex, and potentially reversible, coupling mechanism than previously recognized.

The observed temporal offsets among different pre-seismic parameters can be attributed to the underlying physical mechanisms governing earthquake preparation and energy transfer. Progressive stress accumulation within the lithosphere leads to deformation and instability governed by rock mechanics and fracture processes [73]-[75]. These processes are controlled by fundamental principles of elasticity, thermodynamics, and energy conservation, resulting in the buildup and eventual release of stored strain energy. As a result, different geophysical parameters as thermal, acoustic, and electromagnetic, respond at different times prior to the earthquake, each reflecting specific stages of energy redistribution. This temporal variability highlights the differing sensitivities of these parameters and the complex, multi-layered nature of energy coupling, underscoring the importance of a multi-parametric framework for a comprehensive understanding of seismogenic processes.

Acknowledgements

The authors express their sincere gratitude to Gopi K. Seemala for providing the GPS-TEC software (version 3.5), which was extensively used for the computations in this study. They also acknowledge the NOAA, ECMWF, IGS, GIM, and NASA OMNIWeb data centers for providing the datasets utilized in this work. S. Sasmal acknowledges financial support from the Japan Trust International Research Cooperation Program under NICT for overseas research.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Freund, F. (2011) Pre-Earthquake Signals: Underlying Physical Processes. Journal of Asian Earth Sciences, 41, 383-400.[CrossRef]
[2]	Sasmal, S., Chowdhury, S., Kundu, S., Politis, D.Z., Potirakis, S.M., Balasis, G., et al. (2021) Pre-Seismic Irregularities during the 2020 Samos (Greece) Earthquake (M = 6.9) as Investigated from Multi-Parameter Approach by Ground and Space-Based Techniques. Atmosphere, 12, Article 1059.[CrossRef]
[3]	Pulinets, S. and Ouzounov, D. (2011) Lithosphere-Atmosphere-Ionosphere Coupling (LAIC) Model—An Unified Concept for Earthquake Precursors Validation. Journal of Asian Earth Sciences, 51, 1139-1152.
[4]	Liu, J.Y., Chen, C.H., Lin, C.H., Tsai, H.F., Chen, C.H. and Kamogawa, M. (2011) Ionospheric Disturbances Triggered by the 11 March 2011 M9.0 Tohoku Earthquake. Journal of Geophysical Research: Space Physics, 116, A06319.[CrossRef]
[5]	Hayakawa, M., Molchanov, O.A., Ondoh, T. and Kawai, E. (1996) The Precursory Signature Effect of the Kobe Earthquake on VLF Sub-Ionospheric Signals. Journal of the Communications Research Laboratory, 43, 169-180.
[6]	Clilverd, M.A., Rodger, C.J. and Thomson, N.R. (1999) Investigating Seismoionospheric Effects on a Long Subionospheric Path. Journal of Geophysical Research: Space Physics, 104, 28171-28179.[CrossRef]
[7]	Pulinets, S. and Boyarchuk, K. (2004) Ionospheric Precursors of Earthquakes. Springer.
[8]	Liu, D. (2000) Anomalies Analyses on Satellite Remote Sensing OLR before ChiChi Earthquake of Taiwan Province. Geographic Information Science, 2, 33-36.
[9]	Ghosh, S., Chowdhury, S., Kundu, S., Sasmal, S., Politis, D.Z., Potirakis, S.M., et al. (2022) Unusual Surface Latent Heat Flux Variations and Their Critical Dynamics Revealed before Strong Earthquakes. Entropy, 24, Article 23.[CrossRef] [PubMed]
[10]	Chakraborty, S., Sasmal, S., Chakrabarti, S.K. and Bhattacharya, A. (2018) Observational Signatures of Unusual Outgoing Longwave Radiation (OLR) and Atmospheric Gravity Waves (AGW) as Precursory Effects of May 2015 Nepal Earthquakes. Journal of Geodynamics, 113, 43-51.[CrossRef]
[11]	Gorny, V., Salman, A.G., Tronin, A.A. and Shilin, B.V. (1998) The Earth’s Outgoing IR Radiation as an Indicator of Seismic Activity. Proceedings of the Academy of Sciences of the USSR, 301, 67-69.
[12]	Ghosh, S., Sasmal, S., Maity, S.K., Potirakis, S.M. and Hayakawa, M. (2024) Thermal Anomalies Observed during the Crete Earthquake on 27 September 2021. Geosciences, 14, Article 73.[CrossRef]
[13]	Ouzounov, D. and Freund, F. (2004) Mid-Infrared Emission Prior to Strong Earthquakes Analyzed by Remote Sensing Data. Advances in Space Research, 33, 268-273.[CrossRef]
[14]	Tronin, A.A., Hayakawa, M. and Molchanov, O.A. (2002) Thermal IR Satellite Data Application for Earthquake Research in Japan and China. Journal of Geodynamics, 33, 519-534.[CrossRef]
[15]	Fu, H., Hsiao, L., Yu, X., Lin, C. and Lin, M. (2020) OLR Variability Prior to M ≥ 6 Earthquakes in Taiwan Region. Frontiers in Earth Science, 8, Article 364.
[16]	Saraf, A.K., Choudhury, S., Rawat, V., Das, J.D., et al. (2008) Satellite Detection of Thermal Anomalies Associated with Earthquakes: A Study of Iranian Events. Natural Hazards, 47, 119-135.
[17]	Zhang, H., Zhu, K., Cheng, Y., Marchetti, D., Chen, W., Fan, M., et al. (2023) Atmospheric and Ionospheric Effects of La Palma Volcano 2021 Eruption. Atmosphere, 14, Article 1198.[CrossRef]
[18]	Korepanov, V., Hayakawa, M., Yampolski, Y. and Lizunov, G. (2009) AGW as a Seismo-Ionospheric Coupling Responsible Agent. Physics and Chemistry of the Earth, Parts A/B/C, 34, 485-495.[CrossRef]
[19]	Thurairajah, B., Bailey, S.M., Cullens, C.Y., Hervig, M.E. and Russell, J.M. (2014) Gravity Wave Activity during Recent Stratospheric Sudden Warming Events from SOFIE Temperature Measurements. Journal of Geophysical Research: Atmospheres, 119, 8091-8103.[CrossRef]
[20]	Chakraborty, S., Sasmal, S., Basak, T., Ghosh, S., Palit, S., Chakrabarti, S.K., et al. (2017) Numerical Modeling of Possible Lower Ionospheric Anomalies Associated with Nepal Earthquake in May, 2015. Advances in Space Research, 60, 1787-1796.[CrossRef]
[21]	Yang, S., Asano, T. and Hayakawa, M. (2019) Abnormal Gravity Wave Activity in the Stratosphere Prior to the 2016 Kumamoto Earthquakes. Journal of Geophysical Research: Space Physics, 124, 1410-1425.[CrossRef]
[22]	Piersanti, M., Materassi, M., Battiston, R., Carbone, V., Cicone, A., D’Angelo, G., et al. (2020) Magnetospheric-Ionospheric-Lithospheric Coupling Model. 1: Observations during the 5 August 2018 Bayan Earthquake. Remote Sensing, 12, Article 3299.[CrossRef]
[23]	Hayakawa, M., Schekotov, A., Izutsu, J. and Nickolaenko, A.P. (2019) Seismogenic Effects in ULF/ELF/VLF Electromagnetic Waves. International Journal of Electronics and Applied Research, 6, 1-86.[CrossRef]
[24]	Schekotov, A., Molchanov, O., Hattori, K., Fedorov, E., Gladyshev, V.A., Belyaev, G.G., et al. (2006) Seismo-Ionospheric Depression of the ULF Geomagnetic Fluctuations at Kamchatka and Japan. Physics and Chemistry of the Earth, Parts A/B/C, 31, 313-318.[CrossRef]
[25]	Chowdhury, S., Kundu, S., Ghosh, S., Hayakawa, M., Schekotov, A., Potirakis, S.M., et al. (2022) Direct and Indirect Evidence of Pre-Seismic Electromagnetic Emissions Associated with Two Large Earthquakes in Japan. Natural Hazards, 112, 2403-2432.[CrossRef]
[26]	Hayakawa, M., Schekotov, A., Izutsu, J., Nickolaenko, A.P. and Hobara, Y. (2023) Seismogenic ULF/ELF Wave Phenomena: Recent Advances and Future Perspectives. Open Journal of Earthquake Research, 12, 45-113.[CrossRef]
[27]	Anagnostopoulos, G.C., Karanikola, I. and Marhavilas, P. (2001) Large-Scale Energetic Particle Layers in the High Latitude Jovian Magnetosphere. Planetary and Space Science, 49, 1049-1065.[CrossRef]
[28]	Chakraborty, S., Sasmal, S., Basak, T. and Chakrabarti, S.K. (2019) Comparative Study of Charged Particle Precipitation from Van Allen Radiation Belts as Observed by NOAA Satellites during a Land Earthquake and an Ocean Earthquake. Advances in Space Research, 64, 719-732.[CrossRef]
[29]	Fidani, C. (2015) Particle Precipitation Prior to Large Earthquakes of Both the Sumatra and Philippine Regions: A Statistical Analysis. Journal of Asian Earth Sciences, 114, 384-392.[CrossRef]
[30]	Sasmal, S. and Chakrabarti, S.K. (2009) Ionosperic Anomaly Due to Seismic Activities—Part 1: Calibration of the VLF Signal of VTX 18.2 kHz Station from Kolkata and Deviation during Seismic Events. Natural Hazards and Earth System Sciences, 9, 1403-1408.[CrossRef]
[31]	Kasahara, Y., Muto, F., Horie, T., Yoshida, M., Hayakawa, M., Ohta, K., et al. (2008) On the Statistical Correlation between the Ionospheric Perturbations as Detected by Subionospheric VLF/LF Propagation Anomalies and Earthquakes. Natural Hazards and Earth System Sciences, 8, 653-656.[CrossRef]
[32]	Ray, S., Chakrabarti, S.K., Mondal, S.K. and Sasmal, S. (2011) Ionospheric Anomaly Due to Seismic Activities-III: Correlation between Night Time VLF Amplitude Fluctuations and Effective Magnitudes of Earthquakes in Indian Sub-Continent. Natural Hazards and Earth System Sciences, 11, 2699-2704.[CrossRef]
[33]	Chakrabarti, S.K., Sasmal, S. and Chakrabarti, S. (2010) Ionospheric Anomaly Due to Seismic Activities—Part 2: Evidence from D-Layer Preparation and Disappearance Times. Natural Hazards and Earth System Sciences, 10, 1751-1757.[CrossRef]
[34]	Sasmal, S., Chakrabarti, S.K., Chakrabarti, S. and Chakrabarti, S.K. (2010) Studies of the Correlation between Ionospheric Anomalies and Seismic Activities in the Indian Subcontinent. AIP Conference Proceedings, 1286, 270-290.[CrossRef]
[35]	Maurya, A.K., Venkatesham, K., Tiwari, P., Vijaykumar, K., Singh, R., Singh, A.K., et al. (2016) The 25 April 2015 Nepal Earthquake: Investigation of Precursor in VLF Subionospheric Signal. Journal of Geophysical Research: Space Physics, 121, 10403-10416.[CrossRef]
[36]	Hayakawa, M., Kasahara, Y., Nakamura, T., Hobara, Y., Rozhnoi, A., Solovieva, M., et al. (2010) On the Correlation between Ionospheric Perturbations as Detected by Subionospheric VLF/LF Signals and Earthquakes as Characterized by Seismic Intensity. Journal of Atmospheric and Solar-Terrestrial Physics, 72, 982-987.[CrossRef]
[37]	Sasmal, S., Chakrabarti, S.K. and Ray, S. (2014) Unusual Behavior of Very Low Frequency Signal during the Earthquake at Honshu/Japan on 11 March, 2011. Indian Journal of Physics, 88, 1013-1019.[CrossRef]
[38]	Hayakawa, M., Hobara, Y., Yasuda, Y., Yamaguchi, H., Ohta, K., Izutsu, J., et al. (2012) Possible Precursor to the March 11, 2011, Japan Earthquake: Ionospheric Perturbations as Seen by Subionospheric Very Low Frequency/Low Frequency Propagation. Annals of Geophysics, 55, 95-99.[CrossRef]
[39]	Pal, P., Sasmal, S. and Chakrabarti, S.K. (2017) Studies of Seismo-Ionospheric Correlations Using Anomalies in Phase of Very Low Frequency Signal. Geomatics, Natural Hazards and Risk, 8, 167-176.[CrossRef]
[40]	Politis, D.Z., Potirakis, S.M., Contoyiannis, Y.F., Biswas, S., Sasmal, S. and Hayakawa, M. (2021) Statistical and Criticality Analysis of the Lower Ionosphere Prior to the 30 October 2020 Samos (Greece) Earthquake (M6.9), Based on VLF Electromagnetic Propagation Data as Recorded by a New VLF/LF Receiver Installed in Athens (Greece). Entropy, 23, Article 676. [Google Scholar] [CrossRef] [PubMed]
[41]	Kumar, S., Kumar, S. and Kumar, A. (2022) Earthquakes Associated Subionospheric VLF Anomalies Recorded at Two Low Latitude Stations in the South Pacific Region. Journal of Atmospheric and Solar-Terrestrial Physics, 229, Article 105834.[CrossRef]
[42]	Politis, D.Z., Sasmal, S., Hayakawa, M., Haralambous, H., Datta, A. and Potirakis, S.M. (2024) A Six-Year (2014-2020) Statistical Correlation Study of VLF Terminator Time Shift with Earthquakes in Japan. Remote Sensing, 16, Article 4162.[CrossRef]
[43]	Mannucci, A.J., Wilson, B.D. and Edwards, C.D. (1993) A New Method for Monitoring the Earth’s Ionospheric Total Electron Content Using the GPS Global Network. Proceedings of the Institute of Navigation, Salt Lake City, 20-22 January 1993, 1323-1332.
[44]	Liu, J.Y., Chuo, Y.J., Shan, S.J., Tsai, Y.B., Chen, Y.I., Pulinets, S.A., et al. (2004) Pre-earthquake Ionospheric Anomalies Registered by Continuous GPS TEC Measurements. Annales Geophysicae, 22, 1585-1593.[CrossRef]
[45]	Liu, J.Y., Chen, Y.I., Chuo, Y.J. and Tsai, H.F. (2001) Variations of Ionospheric Total Electron Content during the Chi-Chi Earthquake. Geophysical Research Letters, 28, 1383-1386.[CrossRef]
[46]	Langley, R., Fedrizzi, M., Paula, E., Santos, M. and Komjathy, A. (2002) Mapping the Low Latitude Ionosphere with GPS. GPS World, 13, 41-46.
[47]	Liu, J.Y., Chen, Y.I., Chen, C.H., Liu, C.Y., Chen, C.Y., Nishihashi, M., et al. (2009) Seismoionospheric GPS Total Electron Content Anomalies Observed before the 12 May 2008 Mw7.9 Wenchuan Earthquake. Journal of Geophysical Research: Space Physics, 114, A04320.[CrossRef]
[48]	Ho, Y.Y., Jhuang, H.K., Su, Y.C. and Liu, J.Y. (2013) Seismo-Ionospheric Anomalies in Total Electron Content of the GIM and Electron Density of DEMETER before the 27 February 2010 M8.8 Chile Earthquake. Advances in Space Research, 51, 2309-2315.[CrossRef]
[49]	Kumar, S. and Singh, A.K. (2017) Ionospheric Precursors Observed in TEC Due to Earthquake of Tamenglong on 3 January 2016. Current Science, 113, 795-801.[CrossRef]
[50]	Oikonomou, C., Haralambous, H. and Muslim, B. (2016) Investigation of Ionospheric TEC Precursors Related to the M7.8 Nepal and M8.3 Chile Earthquakes in 2015 Based on Spectral and Statistical Analysis. Natural Hazards, 83, 97-116.[CrossRef]
[51]	De Santis, A., De Franceschi, G., Spogli, L., Perrone, L., Alfonsi, L., Qamili, E., et al. (2015) Geospace Perturbations Induced by the Earth: The State of the Art and Future Trends. Physics and Chemistry of the Earth, Parts A/B/C, 85, 17-33.[CrossRef]
[52]	Sasmal, S., Chowdhury, S., Kundu, S., Ghosh, S., Politis, D., Potirakis, S., et al. (2023) Multi-parametric Study of Seismogenic Anomalies during the 2021 Crete Earthquake (M = 6.0). Annals of Geophysics, 66, SE646. [Google Scholar] [CrossRef]
[53]	Ghosh, S., Sasmal, S., Midya, S.K. and Chakrabarti, S.K. (2017) Unusual Change in Critical Frequency of F2 Layer during and Prior to Earthquakes. Open Journal of Earthquake Research, 6, 191-203.[CrossRef]
[54]	Yang, S.S., Potirakis, S.M., Sasmal, S. and Hayakawa, M. (2020) Natural Time Analysis of Global Navigation Satellite System Surface Deformation: The Case of the 2016 Kumamoto Earthquakes. Entropy, 22, Article 674.[CrossRef] [PubMed]
[55]	Ghosh, S., Sasmal, S., Naja, M., Potirakis, S. and Hayakawa, M. (2022) Study of Aerosol Anomaly Associated with Large Earthquakes. Advances in Space Research, 71, 129-143.
[56]	Hayakawa, M., Izutsu, J., Schekotov, A., Yang, S., Solovieva, M. and Budilova, E. (2021) Lithosphere-Atmosphere-Ionosphere Coupling Effects Based on Multiparameter Precursor Observations for February-March 2021 Earthquakes (M~7) in the Offshore of Tohoku Area of Japan. Geosciences, 11, Article 481.[CrossRef]
[57]	Wu, L., Qi, Y., Mao, W., Lu, J., Ding, Y., Peng, B., et al. (2023) Scrutinizing and Rooting the Multiple Anomalies of Nepal Earthquake Sequence in 2015 with the Deviation-Time-Space Criterion and Homologous Lithosphere-Coversphere-Atmosphere-Ionosphere Coupling Physics. Natural Hazards and Earth System Sciences, 23, 231-249.[CrossRef]
[58]	Marchetti, D., De Santis, A., D’Arcangelo, S., Poggio, F., Piscini, A., A. Campuzano, S., et al. (2019) Pre-Earthquake Chain Processes Detected from Ground to Satellite Altitude in Preparation of the 2016-2017 Seismic Sequence in Central Italy. Remote Sensing of Environment, 229, 93-99.[CrossRef]
[59]	Ouzounov, D., Pulinets, S., Davidenko, D., Rozhnoi, A., Solovieva, M., Fedun, V., et al. (2021) Transient Effects in Atmosphere and Ionosphere Preceding the 2015 M7.8 and M7.3 Gorkha-Nepal Earthquakes. Frontiers in Earth Science, 9, Article 757358.[CrossRef]
[60]	Chetia, T., Sharma, G., Dey, C. and Raju, P.L.N. (2020) Multi-Parametric Approach for Earthquake Precursor Detection in Assam Valley (Eastern Himalaya, India) Using Satellite and Ground Observation Data. Geotectonics, 54, 83-96.[CrossRef]
[61]	Hayakawa, M., Schekotov, A., Izutsu, J., Yang, S., Solovieva, M. and Hobara, Y. (2022) Multi-Parameter Observations of Seismogenic Phenomena Related to the Tokyo Earthquake (M = 5.9) on 7 October 2021. Geosciences, 12, Article 265.[CrossRef]
[62]	Zhang, Y., Wang, T., Chen, W., Zhu, K., Marchetti, D., Cheng, Y., et al. (2023) Are There One or More Geophysical Coupling Mechanisms before Earthquakes? The Case Study of Lushan (China) 2013. Remote Sensing, 15, Article 1521.[CrossRef]
[63]	Dobrovolsky, I.P., Zubkov, S.I. and Miachkin, V.I. (1979) Estimation of the Size of Earthquake Preparation Zones. Pure and Applied Geophysics, 117, 1025-1044.[CrossRef]
[64]	Fritts, D.C. and Alexander, M.J. (2003) Gravity Wave Dynamics and Effects in the Middle Atmosphere. Reviews of Geophysics, 41, 1-64.[CrossRef]
[65]	Seemala, G.K. and Valladares, C.E. (2011) Statistics of Total Electron Content Depletions Observed over the South American Continent for the Year 2008. Radio Science, 46, RS5019.[CrossRef]
[66]	Seemala, G.K. (2023) Estimation of Ionospheric Total Electron Content (TEC) from GNSS Observations. In: Atmospheric Remote Sensing, Elsevier, 63-84.[CrossRef]
[67]	Klotz, S. and Johnson, N.L. (1983) Encyclopedia of Statistical Sciences. Wiley.
[68]	Liu, J., Zhang, X., Novikov, V. and Shen, X. (2016) Variations of Ionospheric Plasma at Different Altitudes before the 2005 Sumatra Indonesia Ms 7.2 Earthquake. Journal of Geophysical Research: Space Physics, 121, 9179-9187.[CrossRef]
[69]	Pulinets, S., Krankowski, A., Hernandez-Pajares, M., Marra, S., Cherniak, I., Zakharenkova, I., et al. (2021) Ionosphere Sounding for Pre-Seismic Anomalies Identification (INSPIRE): Results of the Project and Perspectives for the Short-Term Earthquake Forecast. Frontiers in Earth Science, 9, Article 610193.[CrossRef]
[70]	Tariq, M.A., Shah, M., Hernández-Pajares, M. and Iqbal, T. (2019) Pre-Earthquake Ionospheric Anomalies before Three Major Earthquakes by GPS-TEC and GIM-TEC Data during 2015-2017. Advances in Space Research, 63, 2088-2099.[CrossRef]
[71]	Hayakawa, M., Solovieva, M., Kopylova, G., Hirooka, S., Sasmal, S., Nanda, K., et al. (2025) Multi-Parameter and Multi-Layer Observations of Electromagnetic Precursors to a Huge Hokkaido Earthquake (M = 6.7) on 5 September, 2018, and Lithosphere-Atmosphere-Ionosphere Coupling Channel. Atmosphere, 16, Article 1372.[CrossRef]
[72]	Sasmal, S., Nanda, K., Hayakawa, M., Solovieva, M., Kopylova, G. and Potirakis, S.M. (2025) Towards Understanding Earthquake Preparatory Dynamics: A Multi-Parametric Investigation of the 2025 Kamchatka Mw 8.8 Event. Atmosphere, 16, Article 1328.[CrossRef]
[73]	Scholz, C.H. (2002). The Mechanics of Earthquakes and Faulting. 2nd Edition, Cambridge University Press. [Google Scholar] [CrossRef]
[74]	Brace, W.F. and Byerlee, J.D. (1966) Stick-Slip as a Mechanism for Earthquakes. Science, 153, 990-992.[CrossRef] [PubMed]
[75]	Kanamori, H. (1977) The Energy Release in Great Earthquakes. Journal of Geophysical Research, 82, 2981-2987.[CrossRef]