Numerical Modeling of Coseismic Tropospheric Disturbances Arising from the Unstable Acoustic Gravity Wave Energetics

: Numerous recent studies report the Coseismic Tropospheric Disturbances (CTDs) during large earthquakes. Their presence suggests the importance of atmospheric seismology in a possible earthquake forecasting scenario. The origin mechanism and associated energetics of CTDs are not well understood though the observations associate them with the atmospheric waves. We present the numerical modeling of coupled dynamics of Gravity waves (GWs) and convective instability (CI) in the dry troposphere that produces the CTDs, in the form of pressure disturbances, of observed magnitudes. The study reveals the altitude and epicentral distribution of CTDs and elaborates the relative role of GWs and CI in the generation and intensiﬁcation of CTDs. The study ﬁnds that mega and strong earthquakes disturb the troposphere to a similar level as the severe meteorological weather. However, CTDs become stronger in the conditionally unstable troposphere owing to the convectively unstable wave energetics. The present study also reveals the formation of CTDs at off-epicenter distances which is due to the epicentral inhomogeneity of the ground vibration and subsequent nonlinear feedback that results in the atmospheric inhomogeneity and the excitation of secondary GWs.


Introduction
Troposphere hosts varieties of disturbances during near-surface meteorological, orographical, and lithospheric or seismic activities [1]. Among these activities, the former two forcings are omnipresent, the later forcing limits around the epicenter and depends on the magnitude of the earthquakes. Several studies in recent years report the detections of Coseismic Tropospheric Disturbances (CTDs) during strong earthquakes of magnitudes larger than 7 [2,3]. These CTDs are observed as tropospheric delay in the radio signal, as registered from the GNSS receivers scattered across the epicenter [4,5]. Their presence and associated delay in the radio signal influences the accuracy of the navigation and communication system. They are observed concurrently with the tropospheric pressure disturbances that are distinguished from the nonseismic components based on their association to the epicenter [6].
The origin mechanism and associated energetics of coseismic atmospheric disturbances at the ground are well understood. In the case of atmospheric disturbances of meteorological and orographical origins, Acoustic Gravity waves have received increasing attention, primarily due to their role in the generation of small to mesoscale structures. Mikumo (1968) [7] invoked energetics of Acoustic Gravity Waves (AGWs) to explain the detection of coseismic atmospheric disturbances at the ground during the great Alaskan earthquake of magnitude 9. He showed that most of the observed features of ground pressure disturbances are the manifestation of the AGWs of seismic origin. Based on infrasonic monitoring, numerous studies report the detection of coseismic pressure disturbances that may have their origin at various atmospheric heights, notably at stratospheric and thermospheric reflecting layers [8]. The presence of coseismic pressure disturbances at the long epicentral distance of about 6000 km is owing to the robust energetics of AGWs in the reflecting layers [2].
On the other hand, the detection of CTDs and their association to the energetics of AGWs are relatively recent topics in atmospheric seismology [6]. On the theoretical ground, the origin mechanism and origin altitudes of CTDs are the topics that remain unattended and unresolved to date. Recent observations suggest that the CTDs are the result of Lithosphere-Atmosphere coupling [4,5] that eventually leads to the Coseismic atmospheric Disturbances in the middle and upper atmosphere [9]. The coupling is established and energized by the AGWs in response to the ground uplift during the earthquakes. Owing to the decrease in density (or increase in buoyancy) with the height, the amplitude of AGWs from tiny ground uplift amplifies to many order at the thermospheric heights and forms the detectable atmospheric disturbances [9,10]. In the context of CTDs, the role of robust energetics of the AGWs is not examined. In particular, in the troposphere, where buoyancy increases rapidly with height owing to the steepest density gradient, the amplification may attain a factor of about 10-100 and therefore, CTDs of significant magnitudes are expected to be excited. Moreover, the convectively unstable troposphere is likely to intensify the CTDs. The present theoretical study examines the formation of CTDs from the lithosphere-atmosphere coupling energized by the gravity waves and convective instability. Moreover, the study examines the role of ambient conditions in the intensification and spatial distribution of CTDs.

Materials and Methods
We employ a nonhydrostatic and nonlinear simulation model of AGWs [10,11]. The model solves the following governing equations of wave amplitude (u), atmospheric density (ρ) and pressure (p): Here, Π ν and Π nl are the time derivatives of the viscous and nonlinear forces, ν = µ/ρ and η = η/ρ are the first and second kinematic viscosities, respectively, (µ, η) are the first and second dynamic viscosities respectively and γ is the specific heat ratio.
Much of the insight of AGWs can be obtained from a simplified form of the wave Equation (1), as presented in Appendices A.1 and A.2. The derivation considers strong vertical (y) gradient condition that prevails in the atmosphere and assumes nonlinearnonlocal plane-wave solution of the following form: The derivation also employs the method of separation of variables, i.e., u y = u ys (y)u yt (t), u x = u xs (y)u xt (t) We derive a set of simplified equations for the vertical amplitude of AGWs in the nonviscous atmosphere of the following form: where, Here, Ω b is the nonisothermal nonhydrostatic Brunt-Vaisala frequency (Equation (6.7a) of [12]), γ e and γ ad are the lapse rate and adiabatic lapse rate respectively, γ is the ratio of the specific heats and c is the sound speed. The conditions Ω 2 b > 0 and Ω 2 b < 0 i.e., |γ e | < |γ ad | and |γ e | > |γ ad |, correspond to the stable and unstable conditions respectively. The condition Ω 2 b > 0 excites pure AGWs and the condition Ω 2 b < 0 excites CI. In the troposphere, both conditions may prevail at different altitudes. Therefore, the troposphere may host both AGWs and CI. The present study excludes the compressional dynamics and considers the energetics of gravity waves (GWs) and convective instability. In this case, the pressure disturbance is due to the air motion from advection and convection.
The numerical modeling solves Equations (1)-(3) in the Cartesian coordinate system where x, y and z respectively represents longitude (+ve towards west), altitude and latitude (+ve towards north). The simulation domain covers 0-600 km in altitude, x ep ± 10∆x ep km in longitude and z ep ± ∆z ep km in latitude with identical grid resolution ∆ = 1 km. Here, x ep and z ep are the epicentral coordinates and ∆x ep and ∆z ep determines the size of the fault plane. The simulation begins 30 min before the earthquake onset time (t ep ) and span for 1 h with the the time resolution of 10 s.
Numerous studies have investigated the role of gravity waves in the formation of tropospheric disturbances of nonseismic origins. Three kinds of forcings are responsible for wave energetics: diabatic forcing from condensation/evaporation (e.g., Ref. [13]), obstacle effect associated with the wind shear layer within the troposphere that acts like a mountain and presents effective orographical forcing (e.g., Ref. [14]), and the mechanical oscillator mechanism (e.g., Refs. [15,16]). The present work opts for the mechanical oscillator mechanism in which ground vibration from seismic waves couples mechanically to the atmosphere without the loss of momentum. The continuity of vertical ground velocity across the Earth's surface establishes the coupling. At the lower boundary i.e., at y = 0 km, the continuity of the vertical velocity of the ground associated with the earthquake acts as forcing for the generation of the AGWs, i.e.,  [4] report the detection of CTDs from this earthquake event. Figure 2 demonstrates the ambient atmospheric density and temperature profiles. Figure 3 reveals that two temperature profiles represent convectively stable (|γ e | < |γ ad |) and unstable (|γ e | > |γ ad |) atmospheric conditions. We note that for the convectively unstable temperature profile, the lapse rate is greater than the the adiabatic lapse rate in the altitude range of 2.5-7.5 km.

Results
Based on stable and unstable temperature profiles, we carry out the following numerical experiments:
Unstable+GWs: Simulation of GWs with the unstable temperature profile without the inclusion of CI, 3.
Unstable+GWs+CI: Simulation of GWs with the unstable temperature profile and the inclusion of CI Figure 4 demonstrates the frequency spectrum of the amplitude (u y ) of GWs at an altitude range of 0-20 km. The black circles represent the Brunt-Vaisala period estimated from expression (5) for the stable+GWs case. The spectrum at the ground represents the seismic source spectrum that shows the principal spectral peak at about 2 min and secondary peaks at periods higher than 10 min. We note that the power spectrum at altitudes above ground also reveals spectral peaks at around 2 min and numerous spectral peaks at other periodicities higher than 2 min. Therefore, the GWs sustain the frequency components of the seismic source, as well as generate other frequency components with periods longer than the Brunt-Vaisala period in the troposphere.
The amplitudes u x and u y cause the disturbances in the atmospheric pressure and density, governed by Equations (2) and (3). We define the CTDs of the following form:    Figure 5A, we note the formation of CTDs in the stable troposphere within 20 min from the earthquake onset time. It attains a maximum of about 5 Pascals in the lower troposphere below the altitude of 10 km. At first, long period-long wavefronts with periods longer than 5 min and wavelength of about 10 km are launched in 0-10 km altitudes during the first 25 min from the earthquake onset. Afterward, short period-short wavefronts with periods of about 2 min and wavelengths lower than about 5 km are launched at altitudes higher than 4 km. Therefore, the phase velocities of these wavefronts are less than 40 m/s, a typical range of phase velocity of GWs in the troposphere [16]. The case demonstrates the energetics of GWs to gives rise to the CTDs. Figure 5B reveals the formation of CTDs in the unstable troposphere, similar to the stable+GWs case with few noticeable exceptions. We note the launch of long period-short wavefronts of periods longer than 5 min and wavelength of about 5 km in 0-10 km altitudes during 0-25 min, instead of long period-long wavefronts from the stable+GWs case. The formation of short wavefronts affirms the generation of small scale GWs in the altitude region of the steep unstable gradient. We also note the formation of descending layers in the lower troposphere, a known energetic of GWs. These layers are intense and long-lived in contrast to the weak and short-lived layers in stable+GWs case. In other words, the GWs create coherent descending layers in an unstable troposphere in comparison to incoherent layers in the stable troposphere.  Figure 6 demonstrates the results from the Unstable+GWs+CI numerical experiment. It reveals the formation of descending layers of CTDs, similar to the previous experiments. However, we note the intensification of layers in the altitude range of 2.5-7.5 km where the lapse rate is larger than the adiabatic rate. Therefore, the intensification is owing to the CI that has a positive growth in the unstable altitude range. The condition Ω 2 b < 0 represents the positive growth condition and the growth rate of CI equals to the imaginary part of Ω b . In the case of unstable+GWs without CI, these layers are less resolved and short-lived. Therefore, the unstable+GWs+CI experiment produces the CTDs that persist for a long duration after the mainshock onset time.  Figure 7 demonstrates the evolution of CTDs (pixmaps) and growth rate of CI (contours) in space at a few selected times from the unstable+GWs+CI experiment. The growth rate is positive in the altitude range of 2.5-7.5 km where the ambient lapse rate becomes steeper than the adiabatic lapse rate (Figure 3) The evolution of the growth rate is due to the varying ambient conditions during the production of CTDs. Thus, the evolution reflects the nonlinear dynamics of GWs and CI. We note the development of CTDs at off-epicenter such that they occupy broad epicentral distance, more than the limited ±20 km size of the fault plane. CTDs attain strong amplitudes above the epicenter and in the altitude range of 2.5-7.5 km where the growth rate of CI is positive. The broad epicentral distance coverage of CTDs is owing to the oblique propagating GWs that have their origin from two sources: the epicentral inhomogeneity in the seismic ground forcing, V SISM and secondary generation of GWs from the epicentral pressure gradients developed from the primary GWs.  Figure 8 shows the pressure disturbances at the ground and mean pressure disturbance from the lower troposphere. It reveals the ground pressure disturbance of ∼1 Pa from the ground velocity V SISM ∼0.01 m/s. Moreover, the waveform of the ground disturbance is similar to the waveform of V SISM . Watada et al. (2006) [3] reported the ground pressure disturbances of about 1-5 Pa magnitudes in the vicinity of the epicenter of the Tokachi-Oki earthquake of magnitude 8.3 that had corresponding V SISM ∼0.005 m/s. In addition, he found similar waveforms of ground forcing and ground pressure disturbations for the first 20 min from the earthquake onset. Thus, the observed characteristics of ground pressure disturbances reported by Watada et al. (2006) [3] and simulated characteristics in the present study are similar.

Discussion
We also note that the mean CTDs attain the amplitude of ∼2 Pa in the case of unsta-ble+GWs+CI, in comparison to the lower amplitude in the case of Unstable+GWs case. The contribution from CI begins immediately after the mainshock and continues until 100 min from the mainshock onset. The most significant contribution comes during 40-80 min when the altitude region of the steep gradient is significantly perturbed (Figures 6 and 7), i.e., when the nonlinear dynamics prevail in this altitude region. Therefore, the coupled energetics of GWs and CI can produce intense and long-lasting CTDs. The presence of CTDs alters the refractivity (N) of the troposphere that takes the following form [4]: p T where p and T are in units of millibars and kelvin, respectively. The refractivity introduces the time and phase delays in the radio signal that lead to the following Tropospheric Range Delay or TRD (Λ): where we use the definition of wave velocity v = C/n and the relation between refractivity and refractive index n = 1 + 10 −6 N. Therefore, the residual-refractivity (δN) and residual-TRD (δΛ) can be written as follows: Here, p 0 and T 0 are the pressure and temperature at time t=0. Besides CTDs, the residual-TRD is another variable that can be monitored during the earthquakes [4]. The present study estimates the residual-TRD by estimating residual-refractivity from the pressure and temperature disturbances associated with the CTDs. In the observations, the residual-TRD is estimated from the phase and range delays of dual-frequency radio waves of GNSS satellites [6]. Figure 9 shows the time variation of residual-refractivity and residual-TRD, derived from Equation (7). The altitude range of integration of δN in Equation (7) is 0-20 km. It reveals that similar to the long-lasting CTDs the residual-refractivity and residual-TRD last long after the mainshock. The residual-TRD remains larger than 0.5 mm during the evolution and attains the maximum amplitude of about 2 mm at 20 min from the mainshock onset. Recent studies report the coseismic residual-TRD of 1-3 mm during the Haida Gwaii earthquake of magnitude 7.8 and other similar magnitude earthquakes [4,5]. The residual-TRD of comparable amplitudes from these observations and present simulation affirm that the CTDs of 1-5 Pascals can introduce residual-TRD of a few millimeters in the radio signal. Therefore, the mechanism of unstable+GWs+CI, described in the present study for the generation of CTDs, can explain the observed residual-TRD. Severe and extreme weathers, such as tropical cyclones, produce ground pressure disturbances of about 25 Pa [18]. Mega earthquakes have produced ground pressure disturbances of about 5 Pa [2]. The present numerical study reveals the CTDs of about 5-10 Pa at the tropospheric heights. Therefore, mega and strong earthquakes disturb the troposphere to a similar level as the severe meteorological weather.
The formation of CTDs is a part of atmospheric seismology that reveals their energizations by the GWs and CI. The energizations manifest as the amplification of tiny ground displacement by an order of magnitude at the tropospheric heights and formations of CTDs at a wider epicentral distance than at the ground. Such energizations favor CTDs detection at tropospheric heights during moderate and weak earthquakes that occur before the strong earthquakes. In this regard, the atmospheric seismology of CTDs contributes to the forecasting scenario of strong earthquakes.

Summary
The present study investigates the coupled dynamics of gravity waves and convective instability and their role in the formation of coseismic tropospheric disturbances (CTDs) and Tropospheric Range Delay (TRD). The study presents a number of numerical experiments that employ numerical code of energetics of GWs and CI in response to the ground uplift associated with an earthquake of magnitude ∼8. These experiments reveal the formation of CTDs of magnitude ∼5 Pascals, within about 20 min from the mainshock onset. It is found that the CTDs of 1-5 Pascals can introduce residual-TRD of a few millimeters in the radio signal, similar in magnitude from the recent observations. The CTDs are intense in the region of steep temperature gradients in the lower troposphere. In this altitude region, GWs attain large amplitude and the growth rate of CI is positive. The ground CTDs acquire characteristics similar to the ground pressure disturbances from earlier reports. We argue that the coseismic residual-delay in radio waves from previous reports are owing to the pressure disturbances in the altitude region of steep gradient above the epicenter. The study finds that the GWs alone can produce the CTDs on the conditionally stable troposphere. However, CTDs become stronger in the conditionally unstable troposphere owing to the convectively unstable wave energetics. The present study also reveals the formation of CTDs at off-epicenter distances which is due to the epicentral inhomogeneity of the ground vibration and subsequent nonlinear feedback that results in the atmospheric inhomogeneity and the excitation of secondary GWs.  Data Availability Statement: Simulation data are too large (about 300 MB) in size and we are unable to upload it along with the manuscript. However, if MDPI and/or reviewers demand, we can make them available on googledrive or dropbox.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results. Equations (A5) and (A6) can be solved using the method of separation of variable such that u y = u ys (y)u yt (t), u x = u xs (y)u xt (t) (A7)