A mathematical model of lithosphere–atmosphere coupling for seismic events

Significant evidence of ionosphere disturbance in connection to intense seismic events have been detected since two decades. It is generally believed that the energy transfer can be due to Acoustic Gravity Waves (AGW) excited at ground level by the earthquakes. In spite of the statistical evidence of the detected perturbations, the coupling between lithosphere and atmosphere has not been so far properly explained by an accurate enough model. In this paper, for the first time, we show the result of an analytical-quantitative model that describes how the pressure and density disturbance is generated in the lower atmosphere by the ground motion associated to earthquakes. The direct comparison between observed and modelled vertical profiles of the atmospheric temperature shows the capability of the model to accurately reproduce, with an high statistical significance, the observed temperature fluctuations induced by strong earthquakes.

In the last two decades the observation of ionospheric disturbances that precede and follow earthquakes is one of the most debated topic in the literature. The first detection of changes in ionospheric parameters (e.g., F region critical frequency-foF2) associated to the Alaskan earthquake was reported in 1 . In general, it is well known that the ionosphere is influenced from above by the solar activity, leading to the so-called solar wind-magnetosphere-ionosphere coupling [2][3][4][5] . On the other hand, the ionospheric medium can be influenced from below by atmospheric waves generated in the neutral atmosphere 6 . Since the principal origin of energy in atmosphere is connected to the its most dense (i.e. lowest) layers, it is expected that ionospheric plasma perturbations could be caused by both dynamic tropospheric process, such as cyclones, motion of weather fronts, jet streams 7,8 , and strong tectonic/technogenic sources of atmospheric oscillations, such as earthquakes, tsunami, volcano eruptions etc. [9][10][11] . The hypothesis of causal link between ionospheric perturbation observations associated to seismic activity and neutral atmosphere oscillations, leading to Acoustic Gravity Waves (AGW), was first proposed by 12 . AGW indicates one of the dispersion branches in atmospheric waves characterized by a period of approximately 5 minutes to 10 hours, and wavelength of 10 m to 1000 kilometers. Because of the viscous dissipation of the short wave components, its wavelength increases with altitude, reaching hundreds of meters in the ionospheric D layer (located at an altitude between ∼ 60 and ∼ 90 km) and ∼ 10 km in the F2 layer (located at an altitude of ∼ 400 km). Such waves are the fastest atmospheric oscillations creating upfront of perturbations able to reach the ionosphere 13 . It is worth highlighting that AGWs are able to provide "fast" dynamic coupling between the lower atmosphere and ionosphere, especially those characterized by wavelengths of ∼ 100 − 200 km and phase and group velocities of ∼ 100 − 200 m/s 8 . In addition, during the last two decades there were several publications on experimental data on AGW induced by the earthquake activity (e.g. 14 , and reference therein).
From a theoretical point of view, waves in the atmosphere, as in any other elastic medium, are a relaxing process that appears in response to any perturbation of equilibrium state of the medium. Atmospheric waves may be described using all theoretical provisions of acoustics of fluids and gases. As a consequence, the oscillations in the atmosphere are described by equations that strictly depends on the pressure, inertia forces, and atmospheric vortexes (cyclones and anticyclones), as result of pressure and Coriolis forces, neglecting the inertia term 15,16 .
In this paper we describe a quantitative model that reproduces with high statistical accuracy how pressure and/or density disturbances can be generated in the lower atmosphere by the ground motion associated to an earthquake and how these disturbances can propagate up to the high atmosphere as AGW. As a case study we have www.nature.com/scientificreports/ tested the model using data collected during four earthquakes, making a direct comparison between the observed and modelled vertical profiles of the atmospheric temperature fluctuations induced by the different earthquakes.

Lithosphere: atmosphere coupling model
A seismic event manifests itself through surface waves detected by seismograms, whose dispersion relation is described by Love-Reynolds 17 . The corresponding ground shaking induced by high-magnitude events can excite perturbations at least in the first layer of the atmosphere, that we will call H. Under the hypothesis that the wavelength involved in the perturbation are much grater than H, the dynamics of the upper part of the layer can be roughly described within the shallow water approach 18 . For the sake of simplicity we consider a 1D case, and define η(x, t) as the fluctuation amplitude of the top of the layer and by u(x, t) the horizontal velocity. In the shallow water framework the time evolution of these quantities can be described by two nonlinear equations that can be cast in a conservative form where β is the "batimetry" of the ground (i.e. the impulsive perturbation of the ground), namely the finite shaking of the ground due to the earthquake. The function β(x, t) could be extracted by usual seismograms. However, for the sake of simplicity and without invoking a specific earthquake, we can assume that the impulsive event can be described by a functional shape β(x, t) = β 0 f (x, t)w(t) , where f(x, t) represents the contribution of the seismic surface waves, whose envelope is described by ω(t), related to the finite duration of the earthquake. Then where ω s /k s = v s is the phase speed of Love or Rayleigh surface waves, while, looking at a typical seismograph, we can assume a specific functional shape for the envelope w(t) ∼ t exp(−αt 2 ) , where α −1/2 represents the Strong Motion Duration (SMD) of the seismic event, which is also related to the magnitude.
B y u s i n g t h e a n s at z w h e r e b o t h φ(x, t) = [η(x, t); u(x, t)] a r e p r o p o r t i o n a l t o φ(x, t) = k,ω φ(k, ω) exp[i(kx − ωt)] on the linearized equations (1), and by considering for simplicity that the seismic event generates surface waves by a single component (ω s , k s ) , after some algebra we get the surface perturbation at the first atmospheric layer (see "Excited atmospheric modes during large earthquakes"), which can be cast in a recursive form where the function F is reported in "Excited atmospheric modes during large earthquakes". The solution of (3) allows us to obtain both the normalized amplitude of the perturbation and its dispersion relation (see "Excited atmospheric modes during large earthquakes"). Using Eq. (15) in "Excited atmospheric modes during large earthquakes", we can obtain the wavemodes (k, ω) which can be excited by the earthquake.
Once the fluctuations have been generated roughly at the layer H ⋆ (H* being a characteristic altitude, see "Excited atmospheric modes during large earthquakes"), they give rise to pressure fluctuation. In fact, a parcel of atmosphere at H ⋆ , subject to the vertical displacement η from their equilibrium position, obtained from (15), acquire a vertical velocity w in a way that the Lagrangian pressure fluctuation is given by 19 where p is the Eulerian pressure perturbation at a fixed point in space. By assuming the harmonic ansatz exp[i(k · r − ωt)] for fluctuations, where k = (k x , k y ) , the vertical dependence of both p and w satisfy a set of first-order ordinary differential equations 20 where c 0 is the sound speed, ρ is the mass density, and the intrinsic (Doppler-shifted) wave frequency ω d is defined as ω d = ω − k · u , being u = (u, v) the horizontal velocity. These fluctuations can propagate as an acousticgravity wave through the layered atmosphere, and the waveform at a given height z above H ⋆ can be investigated by using the Wentzel-Kramers-Brillouin (WKB) approximation 21 (see "Propagation of fluctuations through atmosphere: a WKB approach"), thus obtaining the pressure fluctuations at a given height in the atmosphere. As can be seen from Eq. (18) in "Excited atmospheric modes during large earthquakes", the dispersion relation of wave-vectors/frequencies, excited at height, H ⋆ strictly depends on specific earthquake parameters, namely: the Peak Ground Acceleration (PGA), the length of the fault (L), the Strong Motion Duration (SMD -α ), the dominant seismogram frequency ( ω s ) and the phase speed of the surface waves ( v s ).

Discussion
We can investigate the dispersion relation wave-vectors/frequencies excited by a specific earthquake at height H ⋆ .
As case studies we have investigated four earthquakes whose characteristic parameters are reported in Table 1.
We obtained the earthquake parameters from the USGS dedicated website (www.usgs.gov/natural-hazards/ earthquake-hazards/earthquakes). The dispersion relations of fluctuations for each earthquake, in the plane  Figure 1 shows the results obtained for 1995 Kobe earthquake. Looking at the dispersion relation it can be seen that fluctuations η(k, ω) have been roughly excited for wave-vectors ranging from 0.1 ≤ k ≤ 0.35 km −1 , well above k s , and frequencies 0.1 ≤ ω ≤ 0.6 Hz well above ω s . Red dashed line in the panel a) represents the threshold ( ω t = c 0 /h , h being the temperature scale height) for the pressure fluctuations to propagate or do not propagate throughout the atmosphere up to the ionosphere, as purely vertical AGW, i.e. with k · u ≃ 0 . In fact, as a consequence of Eq. (24), waves at frequencies ω excited by the seismic event, which satisfies ω 2 > ω 2 t , are not evanescent and can propagate to high atmosphere. We evaluated ω t = 0.044Hz using the temperature profile retrieved from ERA5, which is the 5 th generation atmospheric data set produced by the European Centre for   11 . The results are presented in panels b), c) and d), displaying vertical atmospheric temperature profile, the atmospheric temperature fluctuations with respect to a 3 km running average and the atmospheric potential energy density, respectively, as obtained from ERA-5 database. A clear AGW propagation, characterized by ∼ 2 km and ∼ 6 km wavelength, is detected in conjunction with the Kobe earthquake occurrence 11 .
Box (B) in Fig. 2 shows the direct comparison between the observed (blue line) and the modelled (red line) temperature fluctuations vertical profile. The latter one can be estimated from the pressure fluctuations model by using the equation of gas in atmosphere p(z) = ρ(z)R * T ′ (z) , and a model profile for the atmospheric density ρ(z) = ρ 0 exp(−γ z) , where R * = R/M , R is the gas constant, M = 28.97 g is the atmospheric mean molecular mass, ρ 0 = 1km/m 3 and γ = 0.06 was atmospheric density mass decay index 23 . The evaluation of the model temperature was done by solving the Eq. (26) in "Propagation of fluctuations through atmosphere: a WKB approach" for wave-vectors and frequencies obtained by the dispersion relation ( Fig. 2 Box A). In particular we used the three pairs corresponding to the maximum values of η/η 0 (see Table 2) and T'(0)=0 as a boundary conditions.
The model is able to correctly reproduce the temperature fluctuation observations with a RMSE (root mean square error) of 0.8 K and a correlation coefficient of 0.86. In addition, to determine whether there is a statistically significant difference between the temperature fluctuations measured and the modelled one, we have performed a χ 2 test (see Table 3 in "χ 2 test"), obtaining χ 2 = 47.3 , suggesting that our model is able to reproduce the observations with > 90% probability. It is interesting to highlight that any phase shift due to a k · u transverse wave-vector increases the χ 2 value. As a consequence, the pure vertical propagation represents the minimum χ 2 condition.
The same analysis has been repeated for the 23 June 2001 Peruvian earthquake. Looking at the dispersion relation (Fig. 3a), we can see that fluctuations η(k, ω) have been roughly excited with wavevectors ranging from 0.1 ≤ k ≤ 0.72 km −1 , well above k s , and frequencies 0.05 ≤ ω ≤ 0.8 Hz well above ω s . Also in this case our model expected an AGW injection in conjunction with the earthquake occurrence, since all the excited modes of η are greater than ω t = 0.051 Hz . This prevision is confirmed by the experimental vertical temperature profile and potential energy density (panels b), c) and d). As for the Kobe event, we modelled the vertical behaviour of the www.nature.com/scientificreports/ temperature fluctuation using the three pairs corresponding to the maximum values of η/η 0 (see Table 2). Even in this case, the model is able to correctly reproduce the temperature fluctuation observations with a RMSE of 0.5 K and a correlation coefficient of 0.88. The χ 2 test confirms again that the model is able to reproduce the observations with > 90% probability with a χ 2 =48.8. Also in this case, we obtained that the pure vertical propagation represents the minimum χ 2 condition. The third event analyzed is the 6 April 2009 L' Aquila earthquake, where ionospheric disturbances were not observed. The analysis of experimental vertical temperature profile (Fig. 4b-d) confirms the lack of AGW injection over the earthquake epicenter. Figure 4a displays the dispersion relation, where η(k, ω) have been roughly excited for wavevectors ranging from 0.28 ≤ k ≤ 0.8 km −1 , well above k s , and frequencies 0.08 ≤ ω ≤ 0.44 Hz well above ω s . However, in this case, all the excited modes are lower than ω t = 0.43Hz : in this case possibly generated AGW are evanescent thus preventing a purely vertical propagation into the high atmosphere.
The analysis of the 19 august 2018 Fiji earthquake shows some peculiar features. In fact, the hypocenter of the seismic event is situated at a large depth, about 300 km, and we observe that fluctuations at the first layer in atmosphere (Fig. 5) have been excited for values of wave-vectors and frequencies lower than the previous cases we analyzed. In particular, the modes excited have wave-vectors of the order of few tens of meters, and frequencies of the order of a few tens of mHz, of the order of the surface waveform k s and ω s (Fig. 5a). The evaluation of ω t ≃ 0.072 Hz, suggests that, in this case, the purely vertical propagation of AGW cannot be excited as a consequence of the earthquake occurrence. However, experimental observations ( Fig. 5b-d) show an AGW propagation concurrently with the earthquake occurrence. In this case the discrepancy with our model prediction is simply related to the absence of a purely vertical propagation. In fact, when we allow a phase shift due to a weak transverse wave-vector, that is k · u ≃ 0.49 Hz, the excited fluctuation (corresponding to the pair: k 1 =0.72 km −1 and ω 1 =0.43 Hz) can propagate as a quasi-vertical AGW, and the temperature profile is almost correctly reproduced by our model even in this particular case (Fig. 5c). In this conditions the RMSE is 2.8K and the correlation coefficient ρ = 0.72 . Finally, for the Fiji earthquake the χ 2 test gives worse result with reference to previous events. In fact, we obtained χ 2 =62.3, corresponding to 65% probability of the model to reproduce the Fiji temperature fluctuations.

Conclusion
The possibility to correctly model, for the first time with high statistical significance, the lithosphere-atmosphere coupling, i.e. how the pressure and density disturbance propagates in the upper atmosphere in case of large earthquakes, is a breakthrough in the co-and post-earthquakes analysis. In particular, modelling how strong www.nature.com/scientificreports/ earthquakes can excite perturbations in the lower atmosphere that propagates vertically as AGW, allows a robust statistical analysis that, in turn, provide a unprecedented tool to trim key parameters to fully reproduce the observed perturbations. In addition, the model allows the a-priori simulation of earthquake effects in the high atmosphere just during and after the event, and opens new scenarios to investigate, through a realistic physical model, co-seismic ionosphere disturbances due to high-magnitude earthquakes. This is a first step towards forecasting earthquake effects in area known to be subjected to large release of energy triggered by tectonic, volcanic etc, events. More specifically, we have modelled the atmospheric fluctuations excited by a generic seismic event on the top of the first layer of the atmosphere, and estimate its dispersion relation as a function of the characteristic parameters of the earthquake. Then, using the Wentzel-Kramers-Brillouin (WKB) approach 21 , we model the pressure fluctuations of the AGW excited by these near-ground fluctuations. The proposed model is able to provide correct descriptions of AGW emission for three superficial earthquakes. For two of them, where the AGW are expected to propagate and are actually observed 11 , the model provide a direct comparison between computed and observed atmospheric vertical temperature profiles. It shows high reliability, with 90% of probability, to correctly reproduce the temperature fluctuations induced by the earthquake with a low RMSE and correlation coefficient larger than 0.86. The model is also able to correctly reproduce the case when ionosphere disturbances were not observed, i.e. L' Aquila earthquake, showing that AGW should be evanescent.
The choice of an analytical model, despite the possible lower accuracy with respect to a numerical one, allows to control and understand the physical process behind the atmosphere-lithosphere coupling system during active seismic conditions. Finally, as shown in Fig. 6 for the KOBE case, the good statistical agreement between the model and the data will allow to trim the model taking into account different earthquakes features, by the use of a minimum χ 2 technique, while available a large data set, with small statistical error and low systematic uncertainties.

Methods
Excited atmospheric modes during large earthquakes. Let us consider the linearized equations (1) www.nature.com/scientificreports/ where, as well known, it is evident that the batimetry modifies the usual nondispersive oscillations described by the classical relations dispersion ω 2 /k 2 = v 2 0 = Hg . By assuming a flat soil, let us suppose now that a seismic event modify impulsively the "batimetry" β(x, t) , which represents here the shacking of the ground during the earthquake. Introducing the Fourier series where φ means each of the variables η , u and β , from the linearized equations we get a pair of algebraic equations  www.nature.com/scientificreports/ The Fourier coefficient of β(x, t) can be easily calculated so that the sum over the pair (q, r) can be eliminated by the delta functions, and, after anti-transforming, we obtain the two following equations for the Fourier coefficients (10) β(k, ω) = β 0 e i(kx−ωt) π 1/2 4α 3/2 ωe −ω 2 /4α δ k,k s δ ω,ω s Table 3. Temperature fluctuations observed vs modelled for the Kobe earthquake. www.nature.com/scientificreports/ By inserting u(k, ω) into the second equation, we obtain a single equation for the perturbation coefficients which can be cast in the recursive form where the multiplicative factor is defined as By interpreting k and ω in terms of the ground values κ = k/k s and ζ = ω/ω s , we can rewrite eq. (11) as where γ = v s /v 0 , the function ψ contains only parameters of the earthquake www.nature.com/scientificreports/ and β = β 0 /g √ α is a dimensionless parameter. The factorized form of Eq. (13) is particularly useful, because it can be solved recursively starting from a reference value η 0 = η(0, 0) which represents the unperturbed background at the surface eight H. Let us consider a length scale L for wavevectors, and a time scale T for frequencies, then they can be enumerated by using k = (2π/L)n and ω = (2π/T)m , where the pair (n, m) = 0, ±1, ±2, . . . are the integers of the wavevector-frequency plane. The shallow water approach requires L >> H , say kH << 1 . If the earthquake generates surface waves with a length scale , the above condition becomes k/k s << /H . This implies that, for a given lengthscale of surface waves, the perturbation can be generated not so far from the ground and wavevectors can be not so high.
If the earthquake produces surface waves with length scale and time scale P, we have κ n = ( /L)n and ζ m = (P/T)m . Iterating relation (13), the Fourier coefficients at a given wavelength-frequency (κ n , ζ m ) , are then described by where, by fixing a pair (i, j) on the plane and defining we have Dispersion relation. First of all note that the parameter β can be interpreted as the Peak Ground Acceleration (PGA) of the earthquake which, in general, is estimated from seismograms in terms of g. This implies that the function ψ contains the information on the characteristic of earthquake, the only unknown quantity is the height H of the layer of atmosphere we are considering. Since the relation (13) is a recursive equation, we must avoid divergences as well as very low values of the perturbation. We can estimate a characteristic height H ⋆ where the value of the amplitude can induce an efficient perturbation for each mode (k, ω) from the relation ψ(H ⋆ , β, α, v s ) ≃ 1 , that is We can interpret this as follows. If we assume that the earthquake is able to induce perturbations leading to acoustic-gravity waves in the atmosphere, the typical amplitude of the perturbation can be roughly excited at least at a characteristic height of the order of H ⋆ , according to relation (17). Assuming H = H ⋆ , the function Ŵ defined by (16), represents a kind of dispersion relation, in the plane (κ n , ζ m ) , which gives information on the waveforms excited by the seismic event.
As a rough estimate we can use P ∼ ω −1 s , while the time base can be estimate as the earthquake duration, say T ∼ 1/ √ α . On the other hand, since >> H ⋆ , we can assume = rH ⋆ , where r is a free parameter, and L, the length base, can be taken to be as an estimate of the seismic faulting. In this way and the dispersion relation depends on the free parameters r and L. In our calculations of the dispersion relation from real earthquakes, we used r = 10.
Propagation of fluctuations through atmosphere: a WKB approach. In the WKB framework, assuming that the scale height of atmosphere h, the velocity u and the sound speed c 0 are smooth function of z/� , where represents some scale of variation of the variables, as a first-order approximation the Fourier coefficient of pressure fluctuations at a given height z is related to that at H  www.nature.com/scientificreports/ the vertical wavevector is defined as: and the coefficient When the condition of propagation m 2 > 0 is satisfied, either the upper and lower sign in (20) propagate as oblique local plane waves. As m 2 < 0 the wave is evanescent, while m 2 (z) = 0 defines the so-called turning points z where WKB approximation fails. Moreover, WKB approximation diverges at points where either ω d = 0 or ω 2 d = kg . However, neither of these conditions exist for purely vertical propagation when k = 0 21 . If we consider quasi-vertically propagating waves, by assuming k << 1 , from (22) we obtain Purely vertical propagation requires that frequencies ω excited by the seismic event must satisfies ω 2 > c 2 0 /4h 2 , in order to be not evanescent and can propagate to high atmosphere. Under this condition, the pressure fluctuations associated to quasi-vertically propagating waves, according to WKB, are where Equations (25) and (26) represents vertically propagating acoustic-gravity waves, with vertical wavevector (24), generated from the perturbation of frequency ω excited at height H ⋆ by ground motion due to the large earthquake. χ 2 test. Here we show the table for the evaluation of the χ 2 test for the KOBE earthquake event and the comparison between the modelled and the observed atmospheric temperature fluctuations (Fig. 6).