Magnetospheric–Ionospheric–Lithospheric Coupling Model. 1: Observations during the 5 August 2018 Bayan Earthquake

: The short-term prediction of earthquakes is an essential issue connected with human life protection and related social and economic matters. Recent papers have provided some evidence of the link between the lithosphere, lower atmosphere, and ionosphere, even though with marginal statistical evidence. The basic coupling is hypothesized as being via the atmospheric gravity wave (AGW)/acoustic wave (AW) channel. In this paper we analyze a scenario of the low latitude earthquake (Mw = 6.9) which occurred in Indonesia on 5 August 2018, through a multi-instrumental approach, using ground and satellites high quality data. As a result, we derive a new analytical lithospheric–atmospheric–ionospheric–magnetospheric coupling model with the aim to provide quantitative indicators to interpret the observations around 6 h before and at the moment of the earthquake occurrence.


Introduction
In the last few years, there has been an increasing interest in the scientific community towards the short-term forecasting of earthquakes (EQs), since many anomalies, statistically correlated with a seismic activity, have been found in the atmosphere and ionosphere, rather than in the lithosphere [1][2][3]. Ionospheric plasma density perturbations, occurring in its lower-side [4], as well as in its upper-side [5], give the most promising results. As a consequence, a few hypotheses have been introduced in order to justify the coupling among lithosphere, atmosphere and ionosphere. The first one, which is based on the chemical transmission, hypothesizes that the atmospheric conductivity can be perturbed by radon outflow close to the earthquake epicenter (EE), leading to a modification of the atmospheric electric field that drives a variation in the ionospheric plasma density profile [6][7][8]. The second one is based on the emission of atmospheric acoustic gravity waves (AGWs). Namely, such oscillations, developing around the EE, perturb the atmosphere in terms of changes in temperature, pressure, ground motion, etc. AGWs can propagate upward and, thus, can drive disturbances in the ionosphere [9][10][11][12]. The third one is electrostatic transmission. Basically, an electrostatic effect is produced in the lithosphere and released in the lower atmosphere by "stress-induced positive holes" which alter the ionospheric ionization status [13,14]. Nonetheless, because of the insufficient amount of experimental evidence supporting those theories, there are many other aspects of the lithosphere-atmosphere-ionosphere coupling process that cannot be explained yet [3,15,16]. In this scenario, models based on AGWs emission seems to provide the most promising results to correctly understand the lithosphere-atmosphere-ionosphere coupling processes in concomitance with an EQ. In fact, many studies using sub-ionospheric soundings through very low frequency/extreme low frequency band before and during a seismic event support the AGW hypothesis-e.g., in [2,4,9,11,17]. Anyway, it has to be stressed that all these observations give no information on the lithosphere/atmosphere coupling mechanisms, since only the lower ionosphere is studied. As a consequence, in recent years, many works hunted for a possible link in the lithosphere-atmosphere-ionosphere system during active seismic conditions. The authors in Korepanov et al. [18], investigating surface atmospheric pressure and magnetic field data during meteorological events, concluded that AGWs are a possible candidate for seismo-ionospheric coupling. A similar study by Nakamura et al. [19] for the 2004 Niigata-Chuetsu EQ (M6.8), using wavelet analyses, confirmed that enhanced variations in the period of 10-100 min of the surface atmospheric pressure and magnetic field data (which are in the range of AGW) enhanced the lower ionosphere perturbation. The authors in Endo et al. [20], using the Japanese F-net wide-band seismic network data, found that, before the 2007 Niigata-Chuetsu Oki and 2008 Iwate-Miyagi EQs, the ground motions in the AGW ranges were enhanced when the lower ionosphere was perturbed. Moreover, the analysis of the crustal motion based on the Global Positioning System data during the 2011 Tohoku EQ suggested a direct time-correlation with very low frequency/extreme low frequency sub-ionospheric and ultra-low frequency magnetic field variations-e.g., in [21][22][23][24]. To conclude, the literature mentioned above provides some evidence on the coupling among lithosphere, lower atmosphere, and ionosphere via the AGW channel.
The present paper, has a two fold approach: on the one hand, it shows the detailed investigation of a Mw 6.9 low latitude earthquake, registered in Indonesia on 5 August 2018, at 11:46 UT, based on a multi-instrumental analysis, using both ground and satellite data; on the other hand, it provides a new analytical magnetospheric-ionospheric-lithospheric coupling model, allowing us, for the first time, to quantitatively explain the experimental observations ∼ 6 h before and at the moment of the earthquake occurrence. The manuscript is organized in the following way: Section 2 shows all the data and methods used for the entire analysis; Section 3 presents the analysis during the EQ occurrence; Section 4 displays the analysis ∼ 6 h before the EQ; Section 5 discusses the results obtained and the main features of the new analytical model proposed; in Section 6 the conclusions are summarized.

Data and Methods
In the following section, we introduce the data set used and the techniques applied for the analysis.

Atmospheric Temperature and Acoustic Gravity Waves Evaluation
The atmospheric temperature profiles were retrieved from ERA5, which is the 5 th generation atmospheric data set produced by the European Centre for Medium-Range Weather Forecasts [25]. The model produces global and hourly temperature profiles with high resolution (137 different pressure levels) from near surface up to 0.01 hPa (∼80 km altitude) using observations from satellites, radiosondes, dropsondes, aircraft, and radars [25]. The horizontal resolution is ∼0.28 • in both longitude and latitude.
In this study, as a key indicator to estimate the AGW activity, we used the potential energy (E P ), which directly depends on the atmospheric vertical temperature profiles [26,27], that in the specific event were retrieved from ERA5. To evaluate E P values, we used the approach by Yang et al. [28,29]. The potential energy density [27] is defined as: where g = 9.8 m/s 2 is the gravitational acceleration (constant with altitude), N is the Brunt − Väisälä frequency, and T is the fluctuations with respect to the background temperature T. These three variables are functions of altitude z and were derived from the ERA5 profiles. To retrieve T, we made a 2 km moving average. T is obtained by subtracting T from the original temperature profile.
The term T T 2 , representing a variance, was evaluated within a layer of 2-km thickness, using the following equation: where h max and h min represent the top and bottom altitudes of the layer, respectively.

The Vertical Total Electron Content (vTEC)
TEC is the total electron content for the line of sight, from the ground receiver to the GNSS satellite whose instrumentation, tracked by the non-profit university-governed consortium UNAVCO (https://www.unavco.org/), records broadcast signals from the satellite constellation used in this work. Specifically, we selected 13 International GNSS Service stations whose geographic coordinates are reported in Table 1, together with the receivers' identifiers, the city and the country. The raw data acquired were converted into standard daily RINEX files suitable for processing. Due to the dispersive nature of the ionosphere it was possible to calculate the TEC starting from the code or carrier phase measurements provided by GNSS receivers on different frequencies. TEC was assumed as the total number of free electrons in a cylinder with cross section of 1 m 2 and height equal to the slant signal path. It was defined as the integral of the electron density (N e ) along the GNSS signal ray path (z), according to the following equation: This quantity can be calculated along each slant path according to the following equation: Equation (4) describes the so-called geometry-free combination in which f 1 and f 2 are the two carrier frequencies of the transmitted signal. In more detail, for GNSS signals f 1 corresponds to L 1 (1575.42 MHz) and f 2 corresponds to L 2 (1227.60 MHz). P 1 and P 2 are the corresponding pseudoranges (i.e., the satellite-receiver optical path), while represents all the frequency-dependent biases induced by the receiver, the satellite and the environment in which the measurement is carried out. As each TEC value obtained from Equation (4) is assimilated to the measure of the electron content along the satellite-receiver path, it is called slant TEC (sTEC). In this way, for each satellite and for each period, TEC measurements take place at different elevation angles and are related to different ionospheric sectors. However, according to the Equation (4), such TEC measurement are affected by biases. So, to minimize such biases, we used the GOPI calibration software (http://seemala.blogspot.it/2011/04/ rinex-gps-tecprogram-version-22.html) that provides calibrated TEC values projected to the vertical (vTEC), by assuming the ionosphere as a single, thin, ionized layer located at 350 km [30]. The vertical projection is allowed to obtain TEC values that do not depend on the position of the GNSS receivers at ground. In this paper, the vTEC time series were calculated over each station in the period between 1 August to 10 August 2018, adopting the method described above. vTEC is expressed in units of 10 16 m −2 = 1TECu. All the results are expressed in Universal Time (UT).

China Seismo-Electromagnetic Satellite (CSES) Data
CSES-01, a sun-synchronous satellite, flies at an altitude of ∼507 km [31], has an orbital descending node time at ∼14:00 LT (local time) and a revisiting period of 5 days. CSES-01 hosts nine payloads: a high-precision magnetometer (HPM) [32], a search coil magnetometer (SCM) [33], an electric field detector (EFD) [34], a Langmuir probe (LAP) [35], a plasma analyzer (PAP) [36], a high energetic particle package (HAPP) [37], an high-energy particle detector (HEPD) [38], a GNSS occultation receiver (GOR) [39], and a tri-band beacon (TBB) [40] to detect magnetospheric signals and perturbations induced from space and ground. In this study we used geomagnetic field (from SCM), electric field (from EFD) and electron density (from LAP) data. The EFD data are in the frequency band from direct current to 2 kHz with a sampling frequency of 2.5 kHz. SCM data are in the frequency band from 10 Hz to 2.2 kHz with a sampling frequency of 51 kHz. Both EFD and SCM vectors were analyzed under the geographical coordinate system. LAP data are sampled every 3 s.

Ground Magnetometer and Magnetospheric Field Line Resonance Frequency Estimation
The magnetometers at ground used for the present analysis come from the INTERMAGNET magnetometer array network, which is a consortium of observatories guaranteeing a common standard data release to the scientific community, leading to possible comparison among measurements at different observation points. In our analysis used 1 s resolution data from Gingin (GNG) (λ = 31.36 • S and φ = 115.71 • E, λ and φ being the geographical latitude and longitude, respectively) and Learmonth (LRM) (λ = 22.22 • S and φ = 114.1 • E).
To evaluate the geomagnetic Field line resonance (FLR) frequency, we used the well established gradient method [41][42][43][44], which relies on geomagnetic observations from a pair of magnetometers slightly separated in latitude and approximately aligned along the same magnetic meridian. The method takes advantage of the fact that the field lines connected to the two stations have similar frequency response, both in amplitude and phase, but are slightly shifted in frequency. In particular, the resonance frequency generally decreases with increasing latitude, except near the plasmapause [45,46] or at very low latitudes (λ < 35 • , [47]). The latitudinal separation between magnetometers must be enough to distinguish the resonance frequency of the two field lines, but, at the same time, ensure sufficient coherency between the two signals. Typically, a separation of 1-3 degrees is required, depending on the latitude and on the quality of the signals. The authors in Waters et al. [44] pointed out that the most reliable determination of the resonance frequency f * is provided by the so-called "cross-phase technique", which, by analyzing the H (North-South) component of the geomagnetic field, identifies f * as the frequency where the phase difference between the equatorward (E) and poleward (P) ground magnetometer observatories φ E − φ P maximizes. The f * detected in this way are assumed to be eigen-frequencies of the toroidal oscillations of the given geomagnetic field line [48][49][50].
It is worth mentioning that GNG and LRM magnetometer stations were the closest to the EE observatories sampled at 1 Hz available for the analyzed period. Anyway, as shown in Menk et al. [46], at low latitudes, the results in evaluating the variation in the eigen-frequencies of the toroidal oscillations of a given geomagnetic field line is substantially unchanged unless the distance between the EE and the observatories is either lower than 20 • in latitude (corresponding to ∼2100 km) , or lower than 25 • in longitude (corresponding to ∼2500 km).

Non-Stationary Signal Decomposition and Their Multiscale Statistical Analysis: The Fast Iterative Filtering Algorithm
The study of techniques for the time-frequency analysis and decomposition of signals is a long lasting line of research which has led over the decades to the development of many important algorithms and approaches [51], which are nowadays commonly used in many research fields. When the signal under analysis is non-stationary, standard methods, such as (Short Time) Fourier Transform and Wavelet Transform, were proven to be inadequate to provide detailed time-frequency information, due to their inherent linearity.
Two decades ago, Huang et al. [52] introduced a game-changing method called Empirical Mode Decomposition (EMD). This is an iterative, local and adaptive data-driven method based on a "divide et impera" approach. The idea is simple, but powerful-we first divide the signal s(t) into several simple components, called Intrinsic Mode Functions (IMFs), plus a trend r(t); then, each IMF is analyzed separately in the time-frequency domain via the computation of the instantaneous frequency of each component [51].
In Huang et al. [52], the IMF is informally defined as an oscillatory function that fulfills two properties: the number of zero crossings equals the number of its extrema, plus or minus one; the envelopes connecting its maxima and minima have to be symmetric with respect to the horizontal axis.
The decomposition produced using the EMD algorithm proved to be successful for a wide range of applications [53,54]. Nevertheless, it is hard to study mathematically (checking, for instance, its convergence) due to the repeated usage of envelopes tailored on the specific signal under study.
In the recent years the so-called Iterative Filtering algorithm and its generalizations have been proposed [55][56][57][58], which are based on iterations and do not require any "a priori" assumption on the signal under analysis.
The structure of IF resemble the EMD one. The key difference is in the signal moving average computed in IF, which is obtained as the convolution of the signal s(t) with an a priori chosen filter function w(t). This apparently small difference between the way IF and EMD compute the moving average opened the doors to the mathematical analysis of IF [53,56,57,[59][60][61] such as the demonstration of its a priori convergence [56,57,61] and its acceleration [59,61] in what is called the Fast Iterative Filtering (FIF) method. FIF allows us to compute the exact same decomposition as IF when the signal is extended periodically at the boundaries. We point out that, whenever the signal is not periodical at the boundaries, we can use the idea proposed in Stallone et al. [54] of pre-extending it as needed, and making it periodical at the newly generated boundaries of the extended signal. FIF proved to be at least two orders of magnitude faster than any other iterative algorithm currently available in the literature for the decomposition of nonstationary signals [59].
For all these reasons, in order to evaluate the multiscale properties of a signal s(t), we first use FIF to decompose s(t) into I MF (t), characterized by a peculiar scale of variability [62], so that where r(t) is the residue of the decomposition. The connection between each IMF and the scale of variability of s(t) has been analyzed by using the Flandrin [51] technique. A dataset characterized by an evident scale separation can be decomposed into two contributions, namely s(t) = s 0 (t) + δs(t). Here, s 0 (t) represents the baseline while δs(t) is the variation around the baseline. To identify δs(t), we applied the method proposed by Alberti et al. [63] by defining δs(t) as the reconstruction of a subset s 1 of k < m modes, characterized by a standardized mean (i.e., the mean divided by the standard deviation) SM ≈ 0 and by IMF fluctuating at higher frequency. To distinguish between possible instrumental origin fluctuations and real signals, a multiscale statistical analysis is needed. For the different scales , we consider the statistics of the values I MF (t). This technique, called multiscale statistical analysis, calculates the moments of the probability distribution p(I MF (t))-the second (variance σ( )), the third (skewness Sk( )), and the fourth (kurtosis excess K ex = K( − 3)). In addition, it evaluates the relative energy rel , and the Shannon information entropy I( ), respectively defined as: These parameters measure the variability of the statistics of the signal as functions of the scale considered [64]. In particular, K ex ( ) indicates how the different s are rich in rare fluctuations [65]; rel measures how "energetically strong" the component is in Equation (5). I( ) measures the "degree of randomness" of each I MF (t) component of the signal. In our case the scale corresponds to the peculiar frequency of each IMF of the signal s(t).

The Bayan 5 August 2018 Earthquake, Co-seismic Observations
On 5 August 2018, a Mw 6.9 earthquake struck the island of Lombok, Indonesia. The epicenter was located in North Lombok Regency at geographic coordinates of 8.28 • S 116.4 • E (Bayan). The earthquake occurred at 11:46 UT, at a depth of 31.0 km (INGV-datacatalog) and was caused by a shallow thrust fault near the Flores Back Arc Thrust, at which the Australia and Sunda plates under-thrusts the Indonesian volcanic island arcs (USGS data catalog: https://earthquake.usgs.gov/\earthquakes/ eventpage/us1000g3ub/executive). It is worth noticing that this earthquake was a part of the Lombok 2018 sequence characterized by four strong Mw > 6.3 within 1 month, implying a very high-level of seismic activity in the study area. In the following section, we describe all the observations provided during the seismic event. Figure 1 shows the E P calculation for 5 August 2018 at 12:00 UT at the epicenter location. The vertical profiles of the temperature, the background temperature and the temperature deviation are reported in panels (a), (b) and (c), respectively. Figure 1d,e present the squared term of the Brunt − Väisälä frequency and the obtained potential energy, respectively. It can be easily seen that the temperature reaches its absolutely maximum at ∼18 km (black dotted line), corresponding to the tropopause [66]. A similar behaviour can be found in both N 2 and E P value [67]. In addition, the potential energy is sometimes additionally enhanced around the stratopause when the temperature changes rapidly with altitude.

Acoustic Gravity Waves Observations
In general, AGWs injection is revealed in the temperature deviation (Figure 1c), and its wavelength can be calculated by a full period (if any) variation. In this case, a clear wave can be identified between 5 km and 30 km, whose peaks are at 7.6, 14.1, 23.7 and 27.2 km (red dotted lines), respectively. E P were maximized at the same altitudes, confirming the injection of AGW. As a consequence, we estimated a ∼7 km vertical wavelength associated with the AGW. To further check the wave activity around the EE, we evaluated the horizontal distribution of the E p from 3 to 5 August 2018 ( Figure 2). The altitude of these maps was fixed at 27 km, corresponding to the maximum potential energy values between the tropopause (∼17 km) and the stratopause (∼40 km). On 3 August 2018, a clear AGW activity (panel a) is visible at lower latitudes than the EE, caused by a cold front associated with the arrival of a cyclone passing the northern part of Indonesia. On the contrary, 4 August 2018 presents a relatively calm state around the EE (panel b). Interestingly, a strong increase in E P (panel c) is clearly visible on 5 August 2018 over the EE (black circle), with respect to both of the previous days. The horizontal wavelength of the wave activity (7 • -11 • in longitude) was ∼700-1200 km. Such values are consistent with previous results in Gokhberg and Shalimov [68], Mareev et al. [69] and more recently in Yang et al. [29] who found horizontal AGW wavelength associated to EQ or explosions between 400-1000 km.

Vertical Total Electron Content Observations
To identify possible ionospheric anomalies, we determined the variation in vTEC with respect to its background, as follows: • We identified 10 days of August 2018 characterized by both low solar activity (i.e., −10 nT < Sym-H < 5 nT and AE < 100 nT, Sym-H and AE being the geomagnetic disturbed time index [70] and Auroral Electrojet index [71], respectively) and low seismic activity (i.e., M < 2, M being the EQ magnitude) in an area of 3 • × 3 • lat × lon around the EE; • We decomposed the diurnal vTEC observations using the FIF method, which we briefly recalled in Section 2.5. The interested reader can find more details on this algorithm and its pseudo-code in [59,61] FIF code for Matlab is freely available at www.cicone.com); • We evaluated the 10-day average relative energy spectrogram (¯ rel ) after removing the long term trend; •¯ rel is the vTEC background.  According to Galav et al. [72], this may be due to the pre-reversal enhancement, which is attributed to the variation of the vertical F region drifts at low latitudes.    Despite the not-ideal latitudinal distance between the two observatories, we were able to correctly estimate the FLR eigen-frequency behavior ( f * ) throughout the entire day. Indeed, between 01:00 UT and 04:00 UT, a stable maximum of the cross-phase is at 78 ± 2 mHz, as expected [47]. At 5:43 ± 0:05 UT (red dashed line) there is a clear shift of f * at a lower value (71 ± 2 mHz). f * came back to its typical values at 7:36 ± 0:05 UT. Interestingly, in conjunction with the EQ time (green dashed line), there is a more pronounced downward shifting of f * at 62 ± 2 mHz. As expected, in the local night (i.e., between ∼13:00 UT and ∼20:00 UT), the algorithm was not able to identify the FLR eigen-frequency due to the lower ionospheric conductivity values [46].

The Bayan 5 August 2018 Earthquake, Pre-Seismic Observations
This section shows the atmospheric, ionospheric and ground magnetic observations ∼6 h before the EQ occurrence, corresponding to the period of CSES satellite flying over the EE (see Figure S1 in the Supplementary Material.

Atmospheric Temperature Observations
As in the co-seismic analysis section, we evaluated ( Figure 5) the E P values of AGWs at 6:00 UT at the epicenter location. Figure 5a-c show the vertical temperature profile as obtained from ERA5, T as obtained by a 2 km moving average and the deviation T , as computed by subtracting the background from the original temperature profile, respectively.  Figure 5d,e are the squared term of the Brunt − Väisälä frequency and the resulted AGW potential energy. Additionally, in this case, as expected [67], we found that the temperature value reached its absolutely maximum at the tropopause (∼17 km). T shows clear waves identified between 21 km and 37 km, whose peaks are at 21.3, 24.1, and 27.3 km, respectively. The maximization of E P at this altitude confirms the injection of AGW, whose estimated vertical wavelength is ∼3 km.
Additionally, in this case, to confirm the wave activity beyond the EE, we evaluated the horizontal distribution of the E p from 3 to 5 August 2018 ( Figure 6). The altitude of these maps was fixed at 27 km, corresponding to the maximum potential energy values between the tropopause (∼17 km) and the stratopause (∼39 km). Panels (a) and (b) confirm the arrival of a cold front associated to the a cyclone arriving at the northern part of Indonesia on 3 August 2018, and the relative calm state on 4 August 2018. Interestingly, a strong increase in E P (panel c) is clearly visible on 5 August 2018 over the EE (black circle), with respect to both the previous days. The horizontal wavelength of the wave activity is (∼15 • in longitude) ∼1600 km.

CSES Satellite Ionospheric Observations
In order to investigate electromagnetic signals possibly associated with 5 August 2018 earthquake, we used CSES satellite electric field, magnetic field and plasma density data. For this purpose, following the approach used in Bertello et al. [73] for DEMETER data, we first decomposed each signal using FIF, than we calculated the relative energy spectrogram ( rel ). To detect any EM anomaly possibly related to the EQ, we compared each rel to a background defined in the following way: A cell 3 • × 3 • in latitude-longitude centered over the EE, in which we evaluated the time-frequency average¯ rel , is selected. The mean operation is applied only if the ratio R ( f ) = rel rel = 1 ± 4σ( f ), σ( f ) being the standard deviation of¯ rel evaluated for each frequency scale.

4.
Each frequency scale showing a K ex almost null and correspondingly a relative maximum in the Shannon entropy (I) was not considered in the evaluation of the relative energy, since it can be represented as a Gaussian fluctuation characterized by high "degree of randomness"-i.e., instrumental noise [58].
We define the evaluated¯ rel as the background, since it gives a representation of the magnetospheric and ionospheric electromagnetic activity directly induced by the sun. As a consequence, any distinct signal over 1 ± 4σ can be reasonably identified as anomalous.
In the case of the Bayan earthquake, when CSES flew over the EE, the solar activity was very low (1 nT < Sym-H < 6 nT) and AE < 100 nT. As a consequence, we evaluated the¯ rel relative to I k,1 interval. In addition, the background was calculated for both EFD and SCM data, using data from August 2018 to November 2019. Figure 7A shows the electric field background¯ rel,EFD for the geographic North-South (E x , left panel), East-West (E y , middle panel) and vertical (E z , right panel) components. A clear signature of the ionospheric Schumann resonance at ∼7.9 Hz is visible in all components [74]. The peaks detected at a frequency of around 2 Hz are due to the v s × B electric field, induced by the motion of the satellite with a velocity v s into a magnetic field B. The peak around 1 kHz, in both E x and E z , can be related to the signature of the plasmaspheric hiss [75,76]. Finally, the peak around 270 Hz along E x is a portion of the whistler mode chorus generated around L shell = 5 propagating into the plasmasphere [77,78]. Figure 7B shows the relative energy evaluated for 5 August 2018. As expected, the Schumann resonance peak at ∼7.9 Hz is clearly visible in all components. Interestingly, anomalous peaks at ∼180 Hz along E y and E z (magenta vertical dashed lines), and at ∼630 Hz along E y (red vertical dashed lines) can be easily identified with respect to the background.
The results of the same analysis repeated for the magnetic field measurements are presented in Figure 8 along the geographic North-South (B x , left panel), East-West (B y , middle panel) and vertical (B z , right panel) components. Here the¯ rel,SCM background (box A) shows the signature of the second Schumann ionospheric resonance at ∼ 20 Hz. In addition, a peak is clearly visible at 12 kHz, probably related to the lower-hybrid resonance of the ionospheric F2 layer [79]. Figure 8B shows the relative energy of the SCM on 5 August 2018. The second Schumann resonance is clearly visible in all components and an anomalous peak at ∼ 180 Hz appears along B x and B z (magenta vertical dashed lines). Such oscillation has the same frequency and is perpendicular to that detected in the EFD data. As a consequence, it is an electromagnetic wave. The analysis of the Poynting flux (see Figure S3 in the Supplementary Material), confirms that such an electromagnetic wave propagates upward. Figure 9a shows the CSES LAP data (blue line) for the entire selected semi-orbit ( Figure S1 in the Supplementary Material).As expected for solar quiet conditions, the density presents a local maximum around the magnetic equator and greatest variations at high latitudes [80]. In order to find possible low amplitude variations (with respect to the long term trend value) switching on over the EQ cell, we removed the long term trend (red line) and analyzed the plasma density fluctuations (black line) using FIF. Figure 9b,c show the LAP background relative energy evaluated from August 2018 to November 2019 for solar quiet conditions and 5 August 2018 rel , respectively. Two anomalous density components switched on between 5:39 ± 00:02 UT and 05:50 ± 00:02 UT characterized by T 1 = 67 ± 1 s and T 2 = 111 ± 1 s peculiar periods.   9d presents plasma density anomaly (ρ * ) obtained by the superposition of the two frequency components, T 1 and T 2 . It can be easily seen that the ρ * switches on near the EE, reaching its maximum values at EE and then vanishing.

Discussion
The analysis presented in this paper is based on a multi-instrumental observations in relation to the 5 August 2018 Bayan earthquake. We divided the analysis into two parts, relative to the co-seismic and pre-seismic observations, making a direct comparison among atmospheric oscillations, ionospheric plasma and electric field perturbations, and magnetospheric FLR eigenfrequency variations. First of all, atmospheric temperature data confirm the injection, over the EE, of a clear co-seismic AGW characterized by a ∼7 km vertical wavelength. A concomitant (1 min after the EQ occurrence) vTEC perturbation (vTEC ), characterized by a period of ∼97 s, was identified with respect to its background. Finally, a clear decrease in the magnetospheric FLR f * was found ∼2 min after the EQ.
The same analysis was repeated ∼6 h before the earthquake, corresponding to the time of the CSES-01 satellite flying over the EE. Additionally, in this case, at ∼6:00 UT, we detected an AGW characterized by a ∼3 km vertical wavelength. In addition, a vTEC perturbation (vTEC * ), characterized by a period of 112.3 ± 5 s, switches on at ∼5:15 UT and vanishes at ∼6:35 UT. The CSES-01 satellite detected an EM wave, at ∼180 Hz, propagating upwards between the ∼5:40 UT and the ∼5:46 UT. Interestingly, the LAP on-board the satellite identifies anomalous density variations, ρ * , characterized by ∼70 s and ∼110 s time periods. Lastly, a decrease in the magnetospheric f * was found ∼5:43 UT. The observational scenario presented above can be explained in terms of a Magnetospheric-Ionospheric-Lithospheric Coupling (M.I.L.C.) model, based on three causal steps (see Figure 10):

1.
An AGW is generated around the EE, propagating through the atmosphere; 2.
The AGW interacts mechanically with the ionosphere, creating a local instability in the plasma distribution through a pressure gradient. Such plasma variation put the ionosphere into a "meta-stable" state, giving rise, in the E-layer, to a local non-stationary electric current. This, in turn, generates an electromagnetic (EM) wave.

3.
The interaction of such EM waves with the magnetospheric field causes a change in the eigenfrequency of the field line, whose ionospheric footprint is located over the radial projection of the EE.
This picture is supported by the following mathematical descriptions. Starting from the general equations of compressible, inviscid flow under a gravity field, in absence of external forcing, where ρ, p, v are the atmospheric density, pressure and velocity, c s is the sound speed and γ is the adiabatic index, we can easily evaluate the exact expression of an AGW propagating in a stratified non-isothermal atmosphere in terms of unsteady, non-uniform pressure perturbation p related by an adiabatic condition in a convected frame [81].
where A = γg 2 /c s . The solution of Equation (7) is a mechanical wave (Lamb wave) [82,83]. As expected, the dispersion relation of such an excited waveform mainly depends of the phase speed of the surface waves of earthquake v s . The intensity of the perturbation depends on the height of the atmosphere where it was excited, and on the characteristic of the earthquake, namely the phase speed v s , the frequency of surface waves ω s , the Peak Ground Acceleration and the Strong Motion Duration of the earthquake. In this case, the earthquake had a magnitude of Mw = 6.9 and, according to the seismic surface waves dispersion relation, registered a frequency on the ground of ω s 0.58 Hz (USGS data catalog: https://earthquake.usgs.gov/earthquakes/eventpage/us1000g3ub/executive), thus corresponding to a phase speed of c s 4.5 km/s. The measured Peak Ground Acceleration was 0.8 g, and the Strong Motion Duration was of 20 s. According to our model, this induces an intense atmospheric perturbation up to an length of about H 10 km, and the waveform was excited with frequencies in the range 0.5 ≤ ω ≤ 2 Hz and wavevectors 500 ≤ k ≤ 3000 m. This range values agrees with the AGWs detected using ERA-5 data (see Figures 1 and 5). In addition, such perturbation was able to propagate, at the sound speed, in the vertical direction [84] and reach the ionosphere.
Putting the ∇p as the external forcing into the Magneto-Hydro-Dynamic equation of a uniform ionosphere, located at 100 km with respect to the Earth's surface and whose plasma density decrease as z −2 with the altitude z, a ionospheric plasma density variation is induced [85]. Indeed, if the Magneto-Hydro-Dynamic system is solved for the perturbation of the electric field E, an EM wave, in the frequency range between 10-700 Hz, is generated, propagating from the topside ionosphere. Additionally, in this case, such prevision is in agreement with our observations (see Figures 7 and 8). At low latitudes (i.e., 0 • ≤ λ ≤ 30 • ), the concurring contribution of the EM wave energy and of the plasma density variation produces a change in the local FLR eigen-frequency. In fact, the resonance frequencies of a geomagnetic field line, with both ends fixed in the ionosphere depend on the field line length, and on the magnetic field intensity and plasma mass density ρ along the field line [86][87][88][89][90][91][92]. A crude estimation of an FLR eigen-frequency ( f ) is given by the time of flight approximation [93] as: where V A is the Alfvén speed along the field line, B(s) is the magnetospheric field along the field line, and ρ(r) is the density at geocentric distance r. As a consequence, a change in the field line length and/or the V A , caused a variation in the FLR frequency [94]. The result of the modeled f * variation of a field line footprinted at λ = 10 • , under the assumption of a pure dipolar Earth's magnetic field, associated with a pressure gradient driven by an earthquake of Mw = 6.9, is reported in Figure 11.
A clear decrease in f * is visible in coincidence to the pressure gradient ( ∇p) driven by the AGW injected by the earthquake. A similar situation is clearly visible in Figure 4. The complete mathematical description of the model will be presented in the second paper. Figure 11. Behavior of the FLR eigen-frequency as modeled by an atmospheric pressure gradient caused by a seismic event.
Despite the fact that the event described above fits very well with the model presented here, many issues are still opened and need a careful inspection.
First of all, if the detection of an AGW propagating through the atmosphere is quite simple (see Figure 1), its direct association to an EQ phenomena is not trivial. It is well known that the generation of AGWs during an EQ can be due to different reasons, such as unstable thermal anomalies that originate from the outflow of greenhouse gases into the atmosphere in fracture regions of the Earth's crust, the oscillation of the Earth's crust having a block-like structure, and so on. Such anomalies can be found in the atmospheric temperature, pressure and/or conductivity [95]. In any case, earthquakes are not the one and the only driver of gravity waves propagating through the atmosphere. In fact, AGW are generally induced by weather systems, by synoptic-scale atmospheric systems and circulations, and similar processes [28,96,97], and can affect temperature and wind fields in the higher atmosphere. At middle and low latitudes (i.e., where the EQs are more recurrent), AGWs in general are meteorologically excited by convective activities around the cold fronts [28]. In our case, following the approach of Yang et al. [29], in order to be sure that the AGW detected either ∼6 h before and at the moment of 5 August 2018 EQ occurrence, we tried to exclude all the other possibilities. For this purpose we examined the weather reports (https://www.accuweather.com/) for 5 August 2018 and we found a tropical cyclone passing the Indonesia around the 12:00 UT, but the cold fronts associated with such a cyclone were essentially located in northern Indonesia, over the Brunei region (i.e., far from the EE). In addition, the weather data over the Bayan zone showed that thunderstorms occurred on 7 and 12 August. So, we are confident that both the AGWs detected at 6:00 UT and at 12:00 UT can be associated to the seismic activity.
Second, the correct discrimination between ionospheric plasma density variations induced by internal and external origin sources is crucial. In general, both vTEC and plasma density irregularities are directly driven by the solar activity-e.g., in [98]. In the present case, the visual inspection of both the solar wind (SW) parameters and the geomagnetic indices (i.e., Sym-H and AE) confirms that 5 August 2018 was a super-solar quiet day [63,99]. In fact, solar wind parameters (see Figure S2 in Supplementary Material) confirms the absence of any structure coming from the sun. On the other hand, as shown in the section above, the values of Sym-H (−5 nT < Sym-H < 6 nT) and AE (AE<100 nT), confirms the low geomagnetic activity either at low and at high latitudes. As a consequence, we can establish that both vTEC* and vTEC" variations are not driven by the sun and can be reasonably associated with the earthquake activity. For the same reason, we can infer that both the variation in f * and the EM activity detected by the CSES-01 satellite are linked to the seismic activity through the process described in the presented model.
Last but not least, while it is easy to connect the AGWs observed around 12:00 to the earthquake, the possible source exciting the AGW detected around 06:00 UT needs an accurate discussion. In our opinion, the possible explanation of pre-earthquake AGW can be either due to the very high-levels of seismic activity of the Lombok 2018 sequence characterized by four strong Mw> 6.3 within 1 month, or to the outflow of charged gases (e.g., radon) by the Earth's surface [9,29]. Both hypotheses are reliable in the present case. In fact, on the one hand surface deformations generating variations in atmospheric temperature (hence AGW formation) are very common during high seismic activity-e.g., in [4]. Such a possibility could be well represented by the M.I.L.C. model since the physical mechanism beyond is quite similar to the co-seismic situation [6]. On the other hand, the lack of local surface gas emission data (radon or charged gases in general) do not allow us to exclude the second hypothesis. In any case, this kind of discharge is, in general, coupled with sub-ionospheric very low frequency electromagnetic emission [29]. In the present case, the visual inspection of the very low frequency channel in the EFD on-board CSES did not show any possible anomalous signal, while, as shown in Figures 7 and 8, an electromagnetic signal in the extreme low frequency band was clearly detected.

Conclusions
The correct identification of atmospheric, ionospheric and/or magnetospheric perturbations possibly connected to seismic activity has become a research focus in the last few decades. Indeed, thanks to the large increase in ground sounding stations and satellites, many observations of ionospheric and/or atmospheric disturbances, related to seismic activity, have been reported. This paper presents a multi-instrumental analysis of the 5 August 2018 Lombok (Bayan) EQ. There is supported evidence that during the EQ there was an activation of the lithosphere-atmosphere-ionosphere-magnetosphere chain, which, starting from the fault break, generates an AGW able to mechanically perturb the ionospheric plasma density, driving the generation of both EM waves and magnetospheric FLR eigen-frequency variation. The EQ scenario presented here represents an optimal case event to be studied and understood since it happened during both super solar quiet and fair weather conditions. We explained such observations through a new developed analytical "M.I.L.C." model, whose complete analytical description and details will be provided in a forthcoming paper (i.e., paper 2). Interestingly, the observations of the CSES-01 satellite flying over the EE around 6 h before the EQ, confirms both the presence of EM wave activity, coming from the lower ionosphere, and plasma density variation (ρ * ) consistent with the vTEC anomaly detected (vTEC*). Finally, we have to stress that the simultaneous availability of both ground and LEO satellite data represented a peculiar feature of the analyzed EQ event which allows us to test and validate the robustness of the presented M.I.L.C. model.  Figure S1: 5 August 2018 CSES satellite orbit (red line). Green line is the partial dataset analyzed. Red circle is the earthquake epicenter location. Figure S2: Solar Wind parameters on 5 August 2018. From the top to the bottom: Interplanetary magnetic field; North-South component of the Interplanetary magnetic field; Solar Wind speed; Solar wind plasma pressure; Sym-H index. Figure S3: Poynting Flux as a function of latitude, as determined by EFD and HPM on-board CSES satellite along geographical directions. The vertical dashed line represents the earthquake epicenter location.
Author Contributions: M.P. writing-original draft preparation, formal analysis, and methodology; M.M. methodology and validation; R.B. writing-review and editing, investigation, validation and supervision; V.C. writing-review and editing, and investigation; A.C. methodology; G.D. methodology; P.D. resources and data curation; P.U. writing-review and editing; All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.

Acknowledgments:
The authors wish to thank the three anonymous reviewers for their comments and suggestions. In addition, we thank both A. De Santis and G. Consolini for their useful comments and suggestions. We thank the national institutes that support INTERMAGNET for promoting high standards of the magnetic observatory practice (www.intermagnet.org) used in this paper. The authors kindly acknowledge N. Papitashvili and J. King at the National Space Science Data Center of the Goddard Space Flight Center for the use permission of 1-min OMNI data and the NASA CDAWeb team for making these data available. This work made use of the data from the CSES mission (http://www.leos.ac.cn/), a project funded by China National Space Administration and China Earthquake Administration in collaboration with Italian Space Agency and Istituto Nazionale di Fisica Nucleare. The parameters of the 2018 Bayan earthquake are provided by USGS (https://earthquake.usgs.gov/) and INGV (http://terremoti.ingv.it/) datacatalogs. ERA-5 data are processed and carried out by ECMWF within the Copernicus climate Change service (C3S), and the data can be retreived from https://www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era5. The Global GNSS Network (GGN) is operated by UNAVCO Inc. at the location of the Jet Propulsion Laboratory (JPL) for the National Aeronautics and Space Administration (NASA) with support from NASA under NSF Cooperative Agreement No. EAR-1261833, and the data can be retreived from https://www.unavco.org/. A. Cicone is a member of the Italian "Gruppo Nazionale di Calcolo Scientifico" (GNCS) of the Istituto Nazionale di Alta Matematica "Francesco Severi" (INdAM). He thanks the INdAM for the financial support under the "Progetto Premiale FOE 2014" "Strategic Initiatives for the Environment and Security-SIES". M. Piersanti and A. Cicone thank the Italian Space Agency for the financial support under the contract ASI "LIMADOU scienza" n • 2016-16-H0.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: