## Abstract

Global Navigation Satellite System (GNSS) measured Total Electron Content (TEC) is now widely used to study the near and far-field coseismic ionospheric perturbations (CIP). The generation of near field (~500–600 km surrounding an epicenter) CIP is mainly attributed to the coseismic crustal deformation. The azimuthal distribution of near field CIP may contain information on the seismic/tectonic source characteristics of rupture propagation direction and thrust orientations. However, numerous studies cautioned that before deriving the listed source characteristics based on coseismic TEC signatures, the contribution of non-tectonic forcing mechanisms needs to be examined. These mechanisms which are operative at ionospheric altitudes are classified as the i) orientation between the geomagnetic field and tectonically induced atmospheric wave perturbations ii) orientation between the GNSS satellite line of sight (LOS) geometry and coseismic atmospheric wave perturbations and iii) ambient electron density gradients. So far, the combined effects of these mechanisms have not been quantified. We propose a 3D geometrical model, based on acoustic ray tracing in space and time to estimate the combined effects of non-tectonic forcing mechanisms on the manifestations of GNSS measured near field CIP. Further, this model is tested on earthquakes occurring at different latitudes with a view to quickly quantify the collective effects of these mechanisms. We presume that this simple and direct 3D model would induce and enhance a proper perception among the researchers about the tectonic source characteristics derived based on the corresponding ionospheric manifestations.

## Introduction

Imprints of seismic forcing on the ionosphere during earthquake events are extensively studied using Global Positioning System (GPS) – Total Electron Content (TEC) measurement technique^{1,2,3,4,5,6,7}. GPS is part of the Global Navigation Satellite System (GNSS), being operated by the US. The GPS recorded coseismic ionospheric perturbations (CIP), are considered to be very useful in providing the spatio-temporal manifestations of tectonic forcing in the ionosphere. Tracing back to the source of CIP, the CIP evolution surrounding the epicenter highly depends on the source characteristics in terms of crustal deformation pattern which evolves mainly during the rupture propagation^{5,8}.

The ground based measurement of crustal deformation has a stringent requirement of availability of the sounding probe at or very near to the deformed area e.g. GPS receiver in precise positioning mode. With the advent of space based technology like InSAR, it is easy to acquire this information, provided successive satellite passes coincide over the earthquake occurrence region^{9}. But for offshore earthquakes, InSAR also fails in providing reliable measurements of crustal deformation. In other cases, one may have to rely on the model based estimation constrained by the available ground measurements, for e.g. 11 March 2011 Tohoku earthquake^{10}.

Therefore, an alternative tool needs to be explored that can provide reliable information about the tectonic source characteristics. In recent times, efforts are on to identify the tectonic source characteristics, such as rupture propagation direction, crustal deformation pattern and thrust orientations, based on the azimuthal distribution of near field CIP as recorded by GPS^{2,5,6,7,11,12}. So, is it always possible to characterize the tectonic source using GNSS measured near field CIP? It seems possible, but there lays the difficulty in terms of non-tectonic forcing mechanisms that act upon the CIP evolution at ionospheric altitudes^{5,6,7,12}.

The basic physical mechanism responsible for the generation of CIP is the redistribution of ionospheric plasma by tectonically induced neutral atmospheric wave perturbations. These perturbations are mainly induced by the coseismic crustal deformation. The major effective non-tectonic forcing mechanisms at ionospheric altitudes are the (i) orientation between the ambient geomagnetic field and seismic induced neutral wave perturbations^{12,13,14,15} (ii) orientation between the moving satellite line of sights and the wave perturbations^{6,7,16} and (iii) ambient ionospheric electron density gradient. The effects of the geomagnetic field-neutral wave orientation on near field CIP evolution and subsequent North-South asymmetry has been estimated with a realistic simulation using a case study based on the Van earthquake^{12}. Bagiya *et al*.^{6,7} further investigated the CIP efficiency in identifying the tectonic source characteristics, during the 12 May 2015 Nepal and 13 November 2016 Kaikoura earthquakes, by introducing an elementary satellite geometry factor which estimates the wave phase cancellation effects on CIP amplitudes arising from the varying satellite geometry. Based on these studies, it could be stated that GNSS measured near field CIP azimuthal distribution offers a reasonable alternative opportunity to determine the seismic source characteristics from the ionosphere, provided non-tectonic forcing mechanisms favor the evolution of CIP.

Despite these individual case studies, no attempt has been made so far to map the combined effects of the non-tectonic forcing mechanisms on GNSS measured CIP. We propose, for the first time, a simple 3D geometrical model to estimate the collective non-tectonic effects on the evolution of CIP within ~600 km of radial distance from the ionospheric projection of seismic source. The constraint of ~600 km radial distance is explained later in the text. Firstly, we define non-tectonic forcing mechanism involving the geomagnetic field variations as Geomagnetic field - Neutral wave orientation Factor (GNF) and the one which involves the satellite geometry effects as Satellite Geometry Factor (SGF)^{6}. The effects of ambient electron density variations are estimated as Electron Density Factor (EDF).

We first estimate the GNF and SGF individually, based on the propagation of seismo-acoustic wave perturbations in 3D atmospheric space and time^{14,17}, for a Trial earthquake source near the Main Himalayan Thrust (MHT). Since satellite geometry effects significantly depend on the location of the ground receiver, the SGF is estimated for various GPS stations assumed within and ~600 km surrounding the Trial tectonic source. The EDF is estimated based on the TEC variations provided by the International GNSS Service (IGS) TEC maps. We then multiply the GNF, SGF, and EDF and introduce a factor that estimates the collective non-tectonic effects in 3D space for the Trial tectonic source. This factor is termed as the Non-Tectonic Forcing Mechanisms (NTFM) factor and the 3D computation model is captioned as the NTFM model.

We further assess the performance of the proposed NTFM model during actual earthquake events. We consider here five events: (i) 01 April 2014 Iquique earthquake (ii) 16 September 2015 Illapel earthquake (iii) 25 April 2015 Nepal earthquake (iv) 11 April 2012 Sumatra earthquake (Mw 8.6) and (v) 13 November 2016 New Zealand earthquake. The selection of these earthquakes is made very thoughtfully. No studies so far reported the latitudinal variations of non-tectonic forcing mechanisms. To generalize the proposed model application for earthquakes occurring at various geographical locations, we validate the NTFM model using these selected earthquake events that occurred at different latitudes. In a brief, the present approach is first of its kind in proposing a 3D model that can estimate the combined effects of non-tectonic forcing mechanisms on the evolution of GNSS measured CIP at various latitudes. We believe that this simple model would assist the researchers with a reasonably clear notion about the tectonic source characteristics derived based on the corresponding ionospheric perturbations.

## Methodology

### Modeling of seismo-acoustic rays in 3D space and time

The upper atmosphere responds profoundly to the tectonic energy within 1 to 10 mHz. This frequency domain contains both acoustic as well as gravity waves. The acoustic wave propagation is significantly affected by the temperature of a medium. This could be understood from the Eq. (1)^{18} which estimates the acoustic wave velocity. Accordingly, the seismo-acoustic wave propagation in the overlying atmosphere depends highly on the ambient atmospheric temperature.

Here γ is the ratio of specific heat capacities, R is the universal gas constant, T is the temperature and M is the molecular mass density. The varying temperature and density cause the refraction of seismo-acoustic waves with atmospheric altitudes by changing their velocity. We model the propagation of seismo-acoustic rays (k) in space and time based on the wave refraction phenomenon. For this, the neutral atmospheric temperature and density are obtained from the NRLMSISE-00^{19} model. The acoustic rays are modeled up to the peak ionospheric electron density altitude obtained from the IRI-2016 empirical model^{20}.

### 3D model to quantify the effects of non-tectonic forcing mechanisms on the manifestations of GNSS measured CIP

We, firstly, model the Geomagnetic field - Neutral wave orientation Factor (GNF). This factor is similar to the ionospheric coupling factor (α) by Calais *et al*.^{13}. In order to do so, we extract the values of magnetic inclination and declination corresponding to the geographic location of a tectonic source and nearby region based on the IGRF-12 model^{21}. Using these values, we estimate the geomagnetic field unit vector and denote it as ** b**. The unit vector

**and the seismo-acoustic rays, modeled using the method described in (1), are the main ingredients for estimating the GNF. The equation**

*b*^{15,16}to calculate the GNF as follows.

Where λ, φ, and h depict the geographic longitude, geographic latitude and terrestrial altitude respectively. From Eq. (2), the GNF can vary from −1 to 1 (cosine function).

In the second step, we estimate the non-tectonic forcing effects arising due to the ambient electron density variations. For this, we use the International GNSS Service (IGS) TEC maps which provide a snapshot of the global ionosphere at every 15 min, 1 hr, and 2 hrs^{22,23}. These maps are prepared by the IGS iono working group using the global IGS TEC observations (http://cdaweb.gsfc.nasa.gov). In maps, TEC is estimated as vertical TEC at ionospheric piercing point (IPP) altitude of 400 km and at latitude × longitude grid of 2.5° × 5.0° within the spatial domain of −87.5° to 87.5° latitudes and −180° to 180° longitudes. The satellite elevation cutoff is 10°. We extract the TEC values for our preferred geographical grid at 15 min interval and interpolate it to obtain the variations at finer grid points. The interpolated TEC values are then normalized and denoted as, Electron Density Factor, EDF(λ, φ, h_{IPP}). Since the obtained TEC variations are at fixed IPP altitude of 400 km, the h is presented as h_{IPP}.

In the third step, the effects of moving satellite geometry are computed. The orientation between the satellite line of sight (LOS) and vertically propagating seismo-acoustic rays has significant effects in terms of phase integration of the propagating waves. Inspired by Georges and Hooke, (1970)^{6}, Bagiya *et al*.^{16} invoked an elementary version of satellite geometry factor (SGF) to account for the wave phase cancelation effects during varying GPS satellite geometry. This factor was formulated as,

where r is the satellite LOS wave vector and χ is the satellite zenith angle. We modify the algebraic expression for the SGF and reformulate it as follows:

We consider the absolute values of ** k**(λ, φ, h) •

**(λ, φ, h)in Eq. (4) and thus the SGF can vary from 0 to 1. The requirement for the algebraic transformation of Eq. (3) is discussed further in the text. We estimate the SGF for near field GNSS stations for all possible satellite LOS in view. Finally, we integrate the effects of GNF, SGF, and EDF by using the algebraic multiplication and propose a simplified NTFM factor as follows.**

*r*Here × denotes the algebraic multiplication. The computations described under Eqs. (1, 2, 4, and 5) are essentially named as the 3D NTFM model. Since electron density measurements from each ionospheric altitudes over a specific geographical location are not always available, we use the IGS map derived TEC for computing the NTFM factor at various ionospheric altitudes.

### Coseismic crustal displacement fields based on the static solution

The coseismic crustal deformation in terms of horizontal and vertical displacement fields are computed based on daily position estimates of the permanent GPS sites. The GAMIT/GLOBK software package is used to estimate the daily positions of GPS sites^{24,25}. In the GAMIT processing, a loosely constrained daily solution is computed which further processed with GLOBK to obtain an accurate solution. The analysis has been performed using GPS measured positioning time series over 41 days (20 days prior to the event, event day and 20 days after the event) in each case. The difference in position coordinates on the day before and day after the earthquake provides the respective coseismic displacement fields.

### Coseismic Ionospheric Perturbations (CIP) from GPS-TEC observations

GPS data in RINEX format from permanent GPS sites have been analyzed based on the following formula to estimate the slant Total Electron Content (sTEC),

where *f*_{1}and *f*_{2} are the carrier frequencies (1575.42 MHz and 1227.60 MHz respectively), *L*_{1}and *L*_{2} are the carrier phases and *λ*_{1} and *λ*_{2} are the corresponding wavelengths.

The data sampling interval is 30 s. An elevation mask of 30° is applied to exclude the low elevation satellite measurements. The sTEC is further converted to vertical TEC using the mapping function M as follows^{26},

where *R*_{E} denotes the Earth’s mean radius, *E* is the satellite elevation angle, and *h* is the IPP height. To highlight the CIP, the obtained vertical TEC is applied with a bandpass filter of 1 to 10 mHz which removes long period trends in vertical TEC.

### Implantation of 3D NTFM model for Trial earthquake source

We consider a Trial earthquake source at the northern low latitude (Fig. 1) and compute the effects of non-tectonic forcing mechanisms at various ionospheric altitudes over the epicenter through our proposed 3D model. The Trial earthquake is assumed on 02 March 2010 at 12:00 UT. It should be noted that this day was a geomagnetically quiet day. To start with, we estimate the altitudinal variations of acoustic wave velocity in the atmosphere over the Trial source location and present in Fig. S1 along with the neutral atmospheric temperature profile obtained from the NRLMSISE-00 model. Using the estimated acoustic wave velocity profile, we model the propagation of seismo-acoustic rays in space with time based on the method described in (1). Figure 1 demonstrates the modelled propagation of six different rays launched at zenith angles (angle from the vertical) of ~58°, ~38.8°, ~31.6°, ~28.7°, ~27.1°, and ~26.2° which are the respective threshold (i.e. beyond which rays refract downward) angles for ray propagation at ionospheric altitudes of 120 km, 150 km, 200 km, 250 km, 300 km, and 350 km. The inset shows the variation of the threshold angle and maximum horizontal distance as a function of the atmospheric altitudes. It should be noted that the NTFM model computation is performed using the rays modeled at every 0.1° resolution within the estimated threshold angle at a given specific altitude. The 2D manifestation of rays modeled at every 0.1° for all possible launch angles (till ~58° from the zenith) could be seen in Fig. S3 and S4. Further, at an altitude of 350 km, the seismo-acoustic rays can propagate to the radial distance of ~680 km from the ionospheric projection of the Trial earthquake source thus we restrict the proposed model computation for a distance of ~600 km surrounding the tectonic source.

We now compute the GNF, EDF, and SGF. The GNF and SGF are computed based on the obtained ray traces in 3D space and time. The GNF depends on two parameters (Eq. 2), the geomagnetic field and atmospheric temperature through acoustic wave velocity. The acoustic wave velocity, as mentioned, is a function of ambient temperature that varies with time of the day, season, solar activity, and latitude as well. On the contrary, the geomagnetic field exhibits secular variations but it varies significantly with latitudes. This could be verified from Fig. S2 which shows the latitudinal variations of magnetic inclination (dip) over the globe obtained from the IGRF-12 model. The location of the Trial epicenter is also shown in the figure. According to the variations of magnetic inclination, we classify the region between 0º–~ 40° as low, ~ 40°–~ 60° as mid and ~60° and above as high latitude region^{27,28,29}.

Figure 2 shows the GNF estimated at ionospheric peak density altitude of 350 km over the Trial seismic source. The red star in the figure demonstrates the ionospheric projection of the Trial epicenter. The 2D schematic in Fig. S3 represents the interaction between the seismo-acoustic rays and geomagnetic field for the Trial seismic source. Since the source is considered in the northern hemisphere, the geomagnetic field vectors are inclined down from the horizontal as shown in the figure. The angle between ** k** and

**is ~180° towards the equator (case-I) and thus the absolute values of GNF ≅ 1 (Eq. 2). Near to the Trial source, the angle is more than 90° (case-II) while both the vectors become almost perpendicular further poleward (case-III). The GNF approaches 0 value for perpendicular geometry. We present the absolute values of GNF in Fig. 2. It has to be noted that GNF varies from 0 to 1 which shows respectively unfavorable to favorable conditions for the manifestations of ionospheric perturbations from the geomagnetic field geometry point of view. Therefore, GNF favors the equatorward evolution and propagation of CIP in the Trial case.**

*b*The latitudinal gradient of ambient ionospheric density plays an important role in manifesting the amplitude of ionospheric perturbations. In the case of higher electron density, the perturbations can evolve with higher amplitudes. However, the reduced density results in perturbations with smaller amplitudes (e.g. during nighttime). In order to estimate the EDF, we extract TEC from the IGS map with a latitude-longitude grid of 18°N–31°N and 78°E–92°E on 02 March 2010 at 12:00 UT and subsequently interpolate it to obtain the TEC variations at a finer resolution. The EDF is estimated by normalizing the TEC and therefore the variations are positive and less than 1 (Fig. 2). The figure demonstrates that the electron density gradient is relatively higher south of the Trial epicenter and thus the CIP in the south may evolve with higher amplitudes.

It should be noted that TEC in IGS maps is computed at an ionospheric altitude of 400 km. In view of the unavailability of electron density measurements at specific ionospheric altitudes, the obtained TEC variations are considered to represent the spatial variation of ionization density at other altitudes as well. It is understood that though the ionospheric electron density changes with altitude, the latitudinal density gradient more or less remains the same at all altitudes. Therefore, the density gradient derived in terms of the normalized EDF at a given ionospheric altitude may reflect the gradient for other altitudes also. However, during geomagnetic storms the ionospheric – plasmaspheric density changes significantly thus this assumption requires extra caution^{30,31,32}.

Next, we compute the SGF using Eq. (4). The GPS satellites orbit in near-circular paths at ~22,000 km. On account of such high orbital elevation and very less eccentricity (<0.02), the satellite azimuth-elevation angles (i.e. LOS) vary in the same manner at any geographical location at a given time. However, the orientation of satellite LOS with seismo-acoustic rays vary with respect to the station-source location. Thus, the SGF has to be estimated at each ground station adjoining the source. The orientation of a satellite LOS with vertically propagating seismo-acoustic waves can have two extremes. The schematic in Fig. S4 shows realistic interaction between the vertically propagating seismo-acoustic rays and satellite LOS in 2D for station located at 400 km from the Trial source. If both vectors (LOS and wave) are orientated at angles ≅ 180° then the absolute values of SGF ≅ 1 (Eq. 3) and the perturbation phases integrate the CIP amplitudes to zero (case-II in Fig. S4). The favorable satellite LOS geometries occur at orientations near 90° with respect to the wave vector (case-I in Fig. S4) and the SGF values ≅ 0 on the occasions (Eq. 3).

It is noteworthy that the maximum favorable GNF value is 1 while the SGF favors most at 0. This flip-flop behavior is due to the geometry of the respective vectors involved in the computation of these factors. Since we aim to combine both these factors and to propose a single model estimation of non-tectonic forcing mechanisms surrounding a seismic source, the current mathematical expressions of these two factors may not be adequate. Thus, we transform the algebraic expression for SGF to Eq. (4) which computes variations closer to 1 for favorable satellite geometry.

As mentioned, the SGF has to be computed for individual GNSS stations. We assume three stations located at 2°, 4° and 6° south of the Trial source and compute SGF for all possible satellite LOS from a specific station (Fig. 2). From Fig. 2, for a station situated near to the Trial epicenter (200 km), the GNSS satellite geometry moderately supports the evolution of CIP. The SGF gradually becomes favorable as the station-source distance increases. It should be noted that the SGF computation would manifest in the same fashion for a station located in the north, east, and west. The maximum station-source distance is restricted within 6° (~600 km) as beyond this distance, the seismo-acoustic rays start refracting downward from ionospheric peak density altitude of 350 km. Further, since ionospheric peak density altitude is 350 km in the Trial case, we could consider the station till 6°. However, in the case of lower ionospheric peak density altitude, the station has to be closer than this.

In the final step, we estimate the NTFM factor at an ionospheric altitude of 350 km using Eq. (5) and present in Fig. 2. We segregate the model output in three domains: (i) [1 0.5] as highly favorable to favorable, (ii) [0.5 0.3] as favorable to moderately favorable (iii) [0.3 0] as moderately favorable to poor. The NTFM factor, in Fig. 2, suggests that the evolution of ionospheric perturbations is favored equatorward, especially for GPS stations located at 4° to 6° latitudinal distance from the source. The CIP growth might be hampered towards the north mainly due to the poor GNF and SGF. Further, for a station located near to the source (2°), the CIP evolution towards the equator might be moderate due to the moderate SGF. It should be noted that Fig. 2 presents the modeled NTFM factor along with GNF, EDF, and SGF at a single altitude of 350 km.

We extend this analysis for various atmospheric altitudes, from 120 km to 300 km, and present the modeled 3D NTFM factor in Fig. 3. It should be noted that ionospheric altitudes of ~110–120 km and above can only host the ionospheric perturbations due to the fact that ionization density is very small below these altitudes. Therefore, the NTFM factor computation is started at 120 km altitude. The next atmospheric altitude is considered at 150 km and then at every 50 km, the factor is computed for stations located at 2°, 4°, and 6° latitudes from the Trial source (Fig. 3). In Fig. 3, for a station located close to the source (at 2°) the non-tectonic effects are more favorable at lower ionospheric altitudes than the higher ones. Whilst, for a station located at 4° the NTFM factor favors the evolution of ionospheric perturbations at higher ionospheric altitudes. The favorable NTFM factor gradually weakens as the station distance increases beyond 4° from the source. This could be verified from a station located at 6° latitude.

It is important to note that the NTFM factor and their horizontal extent vary significantly with altitudes. The variations in the horizontal extent are attributed to the atmospheric refraction of seismo-acoustic rays (Figs. 1, S3 and S4). Since the GNF and SGF computations are based on the propagation of seismo-acoustic rays, the horizontal extent of the NTFM model output evidently reflects the wave propagation characteristics at various atmospheric altitudes. Further, the altitudinal variations of NTFM factor (i.e. favorable and vice versa) depend on the manifestations of GNF and SGF. The geomagnetic field inclination and declination do not vary significantly with altitudes. But the propagation angles of seismo-acoustic rays vary considerably. Therefore, GNF computed based on these two parameters (Eq. 2) varies with altitudes. For SGF, seismo-acoustic wave vectors, satellite LOS, and station-source distance are the key parameters. In the Trial case, for a station located at 2°, the favorable satellite line of sights happen to be those with higher elevations (83°–90°) that align suitably with the upward propagating seismo-acoustic rays (parallel to wavefronts) and thus favor the evolution of CIP. For stations located at 4° and 6°, the favorable line of sights occur for elevation range of ~46° to ~90° (Fig. S4) and ~33° to ~90° respectively. The varying alignments of such satellite geometries with vertically propagating seismo-acoustic rays significantly control the altitudinal behavior of NTFM factor through SGF.

### Validation

In this section, we validate our proposed NTFM model based on the near field coseismic ionospheric response to actual earthquake events that occurred at different latitudes. For this, we analyze ionospheric variations during the Mw 8.2 01 April 2014 Iquique (IQ henceforth) earthquake, Mw 8.3 16 September 2015 Illapel (IL henceforth) earthquake, Mw 7.8 25 April 2015 Nepal (NP henceforth) earthquake, Mw 8.6, 11 April 2012 Sumatra (SU henceforth) earthquake, and Mw 7.8 13 November 2016 New Zealand (NZ henceforth) earthquake. It should be noted that all these earthquakes occurred on geomagnetically quiescent days.

The IQ and IL earthquakes triggered near the Chile subduction zone. The seismicity over the Chilean region is defined by the subduction of Nazca plate under the South American plate at the Peru-Chile trench. This subduction zone is well known for hosting the great and giant megathrust earthquakes^{33}. The IQ and IL earthquakes are examples of this. The IQ earthquake ruptured the center part of the northern Chile seismic gap (~18°S–24°S)^{34} with hypocenter at 19.61°S 70.77°W at a depth of ~25 km (https://earthquake.usgs.gov). Extensive studies on source characteristics, constrained by various slip models and data suggest that the seismic rupture propagated south of the epicenter during this offshore event^{34,35,36,37,38}.

From a geomagnetic coordinate point of view, the IQ earthquake could be categorized as a *near-equatorial earthquake* (~−15.16° magnetic inclination). Figure S2 could be referred for more information on this. The geographical map in Fig. 4a shows the location of IQ epicenter and the earthquake fault mechanism (thrust). The figure also contains the static horizontal and vertical displacement fields estimated using nearby GPS geodetic observations. The displacement fields, presented with horizontal and vertical arrows, exhibited profound westward and downward movements at the coast. From the figure, the maximum of westward and downward movements respectively of ~57.08 cm and ~15.91 cm were observed at *atjn* station. Further, GPS stations from the south of the epicenter also observed significant displacements. The slip variations along with the estimated displacement fields corroborate the maximum energy propagation south of the epicenter. Therefore, the coseismic ionospheric response should be more intense in the south during this event.

Figure 4b–e demonstrate the ionospheric variations over the IQ source region. GPS satellites with Pseudo Random Noise (PRN) codes 01, 20, and 23 from stations shown in Fig. 4b were found to be suitable to observe the near field CIP. This could be verified from the IPP tracks of these satellites in Fig. 4b. The satellites adequately cover the ionosphere over and around the IQ earthquake fault region. Figure 4c–e show the temporal evolution of CIP as observed by these PRNs. The CIP are derived as detrended vertical TEC. Each time series is presented with same color as that of the respective satellite track. PRN 23 and 01 captured clear and intense CIP far north and south of the epicenter. However, the satellites recorded relatively feeble CIP closer to the epicenter. Interestingly, the CIP as recorded by PRN 01 from far north stations of *areq*, *atic* and *nzca* were more intense. In contrast to these, PRN 20 recorded more CIP amplitudes south of the epicenter. Due to the fact that most of the seismic energy propagated south of the epicenter, the CIP should evolve stronger in the south compared to those in the north. However, maxima of CIP amplitudes, as recorded by PRN 23 and PRN 01 were observed far north of the epicenter.

We attempt to quantify the captured ionospheric imprints of seismic energy during the IQ event in light of the manifestation of non-tectonic forcing mechanisms surrounding the earthquake source region. We run the proposed 3D NTFM model to compute the GNF, EDF, SGF and their combined effects as NTFM factor and present the model output at IPP altitude of 350 km in Fig. 5. The red star demonstrates the ionospheric projection of the IQ epicenter. From the estimated GNF variations, the geomagnetic field favored the CIP evolution more in the north as compared to the south during the IQ earthquake. Moreover, the region far north was relatively more favorable for the evolution of CIP from the geomagnetic field geometry point of view. The EDF exhibited relatively higher density in the south during the earthquake occurrence time. The modelled SGF variations in Fig. 5 suggest that the GPS satellite geometry from stations closer to the IQ epicenter (200 km and less) obstruct the CIP evolution. In contrast to this, the SGF considerably favors the CIP evolution for stations located between 200–400 km, both north and south of the epicenter. Finally, the combined NTFM factor in Fig. 5 suggests that stations located within 200 km of the IQ epicenter were not adequate to capture the seismic energy imprints during the event. However, distant stations could capture this energy in a more profound manner.

We assume that poor manifestations of NTFM factor for the nearby GNSS stations (<200 km) did not allow the maximum released seismic energy to evolve in the ionosphere. But in the north, CIP amplitudes remained higher due to the favorable non-tectonic forcing mechanisms. Significant CIP amplitudes were observed in *areq-PRN 01* and *areq-PRN 23*. From Fig. 5, the NTFM factor exhibited favorable to moderately favorable domain for the evolution of CIP at the IPP locations of these PRNs. Near the IQ epicenter, comparatively lesser CIP amplitudes were observed by *iqqe-PRN 01*, *iqqe-PRN 20*, and *iqqe-PRN 23* which are in accordance with the poor NTFM factor. However, irrespective of a similar or more favorable non-tectonic forcing mechanisms, CIP amplitudes remained smaller from *jrgn* station than those of *areq*. This could be attributed to the larger distance between the ionospheric projection of the IQ epicenter and the IPPs from *jrgn* station. Thus, our proposed 3D NTFM model, in light of the prevailing non-tectonic forcing mechanisms, effectively explains the ionospheric variations induced by the earthquake nucleating at near-equatorial latitude.

The IL earthquake is another great subduction earthquake that ruptured the offshore region of central Chile with a hypocenter located at 31.57°S 71.67°W at a depth of ~22.3 km (https://earthquake.usgs.gov). The geographical map in Fig. 6a shows the IL earthquake epicenter location, fault mechanism (thrust) and static horizontal and vertical displacement fields estimated using the geodetic observations from nearby GPS sites. From the figure, significant trenchward (westward) horizontal displacement is evident at all GPS sites. The maximum horizontal displacement of ~1.43 m was observed at *pfrj* GPS station. A vertical uplift of ~2.67 cm was recorded at *cnba* GPS site while other GPS sites subsided by few millimeters to centimeters. This suggests that the offshore coseismic slip exhibited significant slip distribution below the *pfrj* and *cnba* GPS stations^{39}. From the estimated displacement fields and reproduced slip distribution, it is obvious that most of the seismic moment was released north of the epicenter during this event.

The IL earthquake triggered at the geomagnetic low latitude region. The magnetic inclination at the epicenter location is ~−32.43°. This could be verified from Fig. S2. Thus, we classify this earthquake as a *low latitude earthquake*. The TEC as observed by PRNs 25, 12, and 24 from GPS stations shown in Fig. 6b were scrutinized into for coseismic ionospheric variations during the IL event. The IPP trajectories of these PRNs at the earthquake occurrence time are also shown in the figure. From Fig. 6c–e, all three PRNs could capture discernible CIP surrounding the epicenter. However, higher CIP amplitudes were observed north of the epicenter.

To quantify the effects of non-tectonic forcing mechanisms on the observed CIP manifestations, we run our 3D NTFM model for the earthquake source located in low latitude region. The model output, estimated at IPP altitude of 350 km, in terms of GNF, EDF, SGF, and NTFM factor are presented in Fig. 7. The GNF behavior demonstrates that the geomagnetic field favored the CIP evolution more in the north during the IL earthquake. Also the ambient electron density gradient was higher towards the north. The modelled SGF variations corroborate the earlier results for Trial and IQ earthquakes and suggest that GPS stations located within ~200–400 km of epicenter are more suitable to observe the seismic energy imprints in the ionosphere. The estimated NTFM factors for *cmpn*, *sant* and *pmqe* stations indicate that CIP evolution in the north was highly favored by the non-tectonic forcing mechanisms. However, the CIP could not develop in the south mainly due to the poor GNF and EDF. More clarity on this could be obtained by examining the NTFM factor manifestation along the respective satellite tracks from *cmpn*, *sant* and *pmqe* GPS stations in Fig. 7. From Figs. 6 and 7, the CIP evolution preferentially follows the seismic energy propagation during the IL event.

The NP earthquake ruptured the Main Himalayan Thrust (MHT) and was the largest earthquake to occur in the central Himalayas since the 1934 Bihar-Nepal earthquake. The earthquake nucleated at 28.23°N 84.73°E and at a depth of ~8.2 km on 25 April 2015 at ~06:11 UT (https://earthquake.usgs.gov). Figure 8a shows the location of the NP epicenter and the fault mechanism (thrust). The arrows demonstrate the static horizontal and vertical displacement fields estimated using nearby GPS geodetic observations. The horizontal movement of **~**1.75 m in the south-southwest and uplift of ~1.17 m were observed at *kkn4* station. The *chlm* GPS station subsided by ~0.54 m. From Fig. 8a, most of the seismic energy during the NP earthquake propagated in the east-southeast of the epicenter^{5,9}.

Based on the magnetic inclination variations (Fig. S2), the NP earthquake could be categorized as *low-**mid latitude* event (magnetic inclination ~44.31). From Fig. 8b–e, GPS satellites with PRN codes 16, 23, and 26 were rightly located over the NP fault region during the earthquake period and could capture significant ionospheric signatures from the GPS stations located surrounding the epicenter. It should be noted that PRNs 16, 23, and 26 observed significant CIP amplitudes southeast of the epicenter. IPP tracks of these PRNs in Fig. 8b and temporal evolution of CIP as recorded by respective PRNs in Fig. 8(c–e), specially CIP time series of *ptna-16*, *bhgp-23*, and *ptna-23*, could be verified for this. However, the CIP that observed north and southwest of the epicenter could not exhibit comparable amplitudes.

The observed coseismic ionospheric signatures during the NP earthquake are then examined for the effects of non-tectonic forcing mechanisms. We estimate the GNF, EDF, SGF, and total effects of all three based on the proposed model and present the output at IPP altitude of 350 km in Fig. 9. The GNF favors the CIP evolution south of the NP epicenter i.e. equatorward and impedes the evolution in the north. The GNF estimated during the NP earthquake corroborated well with that of the Trial source. From the estimated EDF over the NP seismic source, the electron density was higher in the southeast of the epicenter. For SGF estimation, we select GPS stations of *chlm*, *smkt*, and *ptna*. The modelled SGF was favorable for the manifestation of ionospheric signatures from stations located at the epicentral distance of ~200 km and beyond (i.e. *smkt* and *ptna*). However, the SGF was moderate to poor from a station closer than this. The estimated SGF also corroborated well with that of the Trial source. The NTFM model output in Fig. 9 demonstrates that non-tectonic forcing mechanisms favored the CIP evolution south of the epicenter during the NP earthquake. Whilst, the CIP in the north were not supported by these mechanisms. Therefore, for a source located in the northern low latitudes, if seismic energy propagates south of the epicenter (equatorward) then the induced ionospheric perturbations could be well captured by the stations located south of the seismic source. It should be noted that though the NTFM factor was favorable south of the NP epicenter, the CIP preferably evolved in the southeast which substantiate the seismic rupture propagation in the same direction.

We extend a similar analysis during the SU (near-equatorial latitude), and NZ (mid latitude) earthquakes. The SU earthquake was one of the largest strike-slip events recorded ever. The rupture process of this earthquake was reported to be quite complex with a sequence of ruptures propagating westward to the epicenter^{40,41}. The strike-slip faults that ruptured were mainly trending west-northwest to east-southeast and north-northwest to south-southwest. Figure S5 shows the coseismic ionospheric response to this event as observed by PRNs 32, 16, and 06 from various Sumatran GPS Array (SuGAR) stations. It could be noted that PRNs orbiting north and south of the epicenter recorded discernible CIP. However, CIP remained quite feeble nearby the epicenter (e.g. *umlh-6, bthl-16*). The NTFM factor along with GNF, EDF, and SGF during the SU event is shown in Fig. S6. Similar to the IQ earthquake, the NTFM factor during the SU event favored the CIP evolution in the north and south of the epicenter and hampered the same nearby the epicenter. Interestingly, the CIP amplitudes were consistently larger towards the north. This corroborated with the fact that maximum vertical displacement during this event was observed in the north of the epicenter^{4,42}. Despite the complex rupture process during the SU event, the coseismic ionospheric response could reflect the maximum tectonic energy propagation reasonably well.

Bagiya *et al*.^{7} analyzed the NZ earthquake and explained the coseismic ionospheric response based on the proposition of tectonic thrust orientations as distinct sources responsible for peculiar CIP azimuthal propagation surrounding the Kaikoura epicenter. This novel proposition was derived considering the non-tectonic forcing mechanisms in the account. However, they could not propose any mechanism to quantify all three non-tectonic effects collectively. It may be recalled that NZ earthquake was a mid latitude event. We reanalyzed this event using our proposed model and considering the Campbell Coseismic Thrust Zone as the source of CIP^{7}. We found that during the NZ event also the non-tectonic forcing mechanisms support the equatorward evolution of seismic induced ionospheric perturbations. This could be realized from the evolution of CIP over the NZ fault region presented in Fig. S7. The NTFM model output during the NZ earthquake is presented in Fig. S8 to acquire better clarity on this.

In the final step, we evaluate the correlation between the peak-to-peak amplitude variations of CIP and corresponding variations of NTFM factor at IPP altitude of 350 km for each satellite and station pairs shown in Figs. 4, 6, 8, S5 and S7 and present in Fig. 10. The correlation analysis is presented for different latitude conditions. Therefore, the figure contains three sets of events analyzed in this case study i.e. near-equatorial, low, and mid latitudes. The region around the ionospheric projection of respective epicenters is divided into various concentric circles with a 100 km span. The correlation is estimated at each range from a radial distance of 100 km and onward. The ionospheric perturbations corresponding to each event are presented with a specific symbol for easy classification.

It could be noted that the correlation between CIP amplitudes and corresponding NTFM factor withstands fairly well at each range of distance. However, the amplitudes weaken beyond 400 km of radial distance in the present case study. We propose that beyond these threshold distances, direct propagation of epicentral energy may not be adequate to induce significant ionospheric signatures even if the NTFM factor is favorable. This allows us to provide a first estimation of the threshold distance for the CIP evolution from the viewpoint of direct epicentral energy propagation and the impact of non-tectonic forcing around the epicenter. The correlation analysis for the radial distance <100 km was not feasible due to a lack of sufficient CIP observations. Further, during the NZ earthquake event, very few CIP were observed after 400 km so no correlation analysis is presented.

From the present study, it should be noted that for earthquakes occurring at near-equatorial latitudes, the NTFM factor favors the evolution of CIP north and south of the epicenter (>200 km) and not nearby the epicenter. The coseismic ionospheric variations during the IQ and SU earthquakes could be verified for this. Further, for earthquake nucleating at low and mid latitudes, the NTFM factor favors the CIP manifestation equatorward. The coseismic ionospheric responses to the NP, IL and NZ earthquakes substantiate this. For easier and quick review of this study, the response of NTFM model for earthquakes occurring at different latitudes is summarized in Table 1.

## Conclusion

A 3D model to map the combined effects of the non-tectonic forcing mechanisms of geomagnetic field, GNSS satellite geometry, and ambient electron density gradient is proposed for the first time. The 3D NTFM model can compute these effects at various ionospheric altitudes depending on the propagation characteristics of seismo-acoustic rays. The model not only successfully explains the ionospheric manifestations during seismic events occurring at different latitudes but also cautions that any correlation between the seismic source manifestations at the ground and corresponding ionospheric perturbations could be erroneous in the absence of quantifying the effects of these mechanisms. Further, the threshold distance from the viewpoint of direct epicentral energy propagation and impact of non-tectonic forcing around the epicenter is modeled for the first time. It should be noted that the proposed 3D model is specifically designed for the spatial analysis of GPS-TEC derived seismic induced ionospheric perturbations.

The preparation of GUI based online version of the proposed model is in progress to facilitate the researcher to quantify the non-tectonic effects on ionospheric perturbations during earthquakes occurring at various latitudes. This model will be shortly online at iigm.res.in. It is believed that the proposed study will assist the researchers to identify the non-tectonic effects on the GNSS measured CIP without performing any further quantitative analysis.

## Data availability

GPS-TEC data over the Chile and Nepal regions are obtained respectively from http://gps.csn.uchile.cl/ and https://www.unavco.org/data/.

## References

- 1.
Calais, E. & Minster, J. B. GPS detection of ionospheric perturbations following the January 17, 1994, northridge earthquake.

*Geophyiscal Res. Lett.***22**, 1045–1048 (1995). - 2.
Astafyeva, E., Lognonné, P. & Rolland, L. First ionospheric images of the seismic fault slip on the example of the Tohoku-oki earthquake.

*Geophys. Res. Lett.***38**, 1–6 (2011). - 3.
Rolland, L. M.

*et al*. The resonant response of the ionosphere imaged after the 2011 off the Pacific coast of Tohoku Earthquake.*Earth, Planets Sp.***63**, 853–857 (2011). - 4.
Sunil, A. S., Bagiya, M. S., Reddy, C. D., Kumar, M. & Ramesh, D. S. Post-seismic ionospheric response to the 11 April 2012 East Indian Ocean doublet earthquake.

*Earth, Planets Sp.***67**, 37 (2015). - 5.
Sunil, A. S.

*et al*. Dependence of near field co-seismic ionospheric perturbations on surface deformations: a case study based on the April, 25 2015 Gorkha Nepal earthquake.*Adv. Sp. Res*. https://doi.org/10.1016/j.asr.2016.11.041 (2016). - 6.
Bagiya, M. S.

*et al*. Efficiency of coseismic ionospheric perturbations in identifying crustal deformation pattern: Case study based on Mw 7.3 May Nepal 2015 earthquake.*J. Geophys. Res. Sp. Phys*.**122**, (2017). - 7.
Bagiya, M. S., Sunil, P. S., Sunil, A. S. & Ramesh, D. S. Co-seismic contortion and coupled nocturnal ionospheric perturbations during 2016 Kaikoura, Mw 7.8 New Zealand earthquake.

*J. Geophys. Res. Sp. Phys*. 1–11, https://doi.org/10.1002/2017JA024584 (2018). - 8.
Otsuka, Y.

*et al*. GPS detection of total electron content variations over Indonesia and Thailand following the 26 December 2004 earthquake.*Earth, Planets Sp.***58**, 159–165 (2006). - 9.
Sreejith, K. M.

*et al*. Coseismic and early postseismic deformation due to the 25 April 2015, Mw 7.8 Gorkha, Nepal, earthquake from InSAR and GPS measurements.*Geophys. Res. Lett*.**43**(2016). - 10.
Bletery, Q.

*et al*. Journal of Geophysical Research: Solid Earth A detailed source model for the M w 9. 0 Tohoku-Oki earthquake reconciling geodesy, seismology, and tsunami records.*J. Geophys. Res. Solid Earth***119**, 7636–7653 (2014). - 11.
Astafyeva, E. I. & Heki, K. Dependence of wave form of near-field coseismic ionospheric disturbances on focal mechanisms.

*Earth Planets Sp.***61**, 939–943 (2009). - 12.
Rolland, L. M.

*et al*. Discriminating the tectonic and non-tectonic contributions in the ionospheric signature of the 2011, Mw7.1, dip-slip Van earthquake, Eastern Turkey.*Geophys. Res. Lett.***40**, 2518–2522 (2013). - 13.
Calais, E. & Minster, J. B. GPS, earthquakes, the ionosphere, and the Space Shuttle.

*Phys. Earth Planet. Inter.***105**, 167–181 (1998). - 14.
Heki, K. & Ping, J. Directivity and apparent velocity of the coseismic ionospheric disturbances observed with a dense GPS array.

*Earth Planet. Sci. Lett.***236**, 845–855 (2005). - 15.
Rolland, L. M., Lognonné, P. & Munekane, H. Detection and modeling of Rayleigh wave induced patterns in the ionosphere.

*J. Geophys. Res. Sp. Phys.***116**, 1–18 (2011). - 16.
Georges, T. M. & Hooke, W. H. Wave-induced fluctuations in ionospheric electron content - A model indicating some observational biases.

*J. Geophys. Res.***75**, 6295–6308 (1970). - 17.
Dautermann, T., Calais, E. & Mattioli, G. S. Global positioning system detection and energy estimation of the ionospheric wave caused by the 13 july 2003 explosion of the soufriere hills volcano, montserrat.

*J. Geophys. Res. Solid Earth***114**, (2009). - 18.
Kinsler, L. E.,

*and A. R. F. Fundamentals of Acoustics*. (John Wiley, Hoboken, N.J., 1962). - 19.
Picone, J. M., Hedin, A. E., Drob, D. P. & Aikin, A. C. NRLMSISE-00 empirical model of the atmosphere: Statistical comparisons and scientific issues.

*J. Geophys. Res. Sp. Phys.***107**, 1–16 (2002). - 20.
Bilitza, D.

*et al*. International Reference Ionosphere 2016: From ionospheric climate to real-time weather predictions. Sp.*Weather***15**, 418–429 (2017). - 21.
Thébault, E.

*et al*. International geomagnetic reference field: The 12th generation international geomagnetic reference field - The twelfth generation.*Earth, Planets Sp*.**67**, (2015). - 22.
Feltens, J. & Schaer, S. IGS products for the ionosphere. In

*IGS position paper, IGS analysis centers workshop ESOC, Darmstadt, Germany, February 9*–*11*, 1998 (1998). - 23.
Wienia, R. J.

*Use of global ionospheric maps for precise point positioning*. (2008). - 24.
Herring, T. A., King, R. W., Mcclusky, S. C. & Sciences, P.

*Introduction to GAMIT / GLOBK. Mass. Instit. Tech*. 1–50 (2015). - 25.
King, R.W. and Bock, Y.

*Documentation for the GAMIT GPS Analysis Software*. (2002). - 26.
Mannucci, A. J., Wison, B. D. & Edwards, C. D. A new method for onitoring the Earth’s ionospheric total electron content using GPS global networks.

*Proc. ION GPS***93**, 1323–1332 (1993). - 27.
Sastri, J. Equatorial anomaly in F-region—A Review. Indian.

*J. Radio Sp. Phys.***19**, 225–240 (1990). - 28.
Hunsucker, R. D. & Hargreaves, J. K.

*The High-Latitude Ionosphere and Its Effects on Radio Propagation*. Cambridge University Press**85**, (2003). - 29.
Rajaram, G. & Pisharoty, P. R.

*The Earth’s Magnetic Field*. (Oxford and IBH Publishing Co. Pvt. Ltd., New Delhi, 1998). - 30.
Piersanti, M., Cesaroni, C., Spogli, L., Alberti, T. & ScienceDirect Does, T. E. C. react to a sudden impulse as a whole? The 2015 Saint Patrick’ s day storm event.

*Adv. Sp. Res.***60**, 1807–1816 (2017). - 31.
Piersanti, M.

*et al*. Comprehensive Analysis of the Geoeffective Solar Event of 21 June 2015: Effects on the Magnetosphere, Plasmasphere, and Ionosphere Systems.*Sol Phys*. https://doi.org/10.1007/s11207-017-1186-0 (2017). - 32.
Pezzopane, M.

*et al*. On some features characterizing the plasmasphere – magnetosphere – ionosphere system during the geomagnetic storm of 27 May 2017.*Earth, Planets Sp.*https://doi.org/10.1186/s40623-019-1056-0 (2019). - 33.
Lomnitz, C. Major Earthquakes of Chile: A Historical Survey, 1535–1960.

*Seismol. Res. Lett.***75**, 368–378 (2004). - 34.
Geersen, J., Ranero, C. R., Barckhausen, U. & Reichert, C. Subducting seamounts control interplate coupling and seismic rupture in the 2014 Iquique earthquake area.

*Nat. Commun.***6**, 6–11 (2015). - 35.
Lay, T., Yue, H., Brodsky, E. E. & An, C. The 1 April 2014 Iquique, Chile, Mw 8.1 earthquake rupture sequence.

*Geophys. Res. Lett.***7**, 1–8 (2014). - 36.
Duputel, Z.

*et al*. The Iquique earthquake sequence of April 2014: Bayesian modeling accounting for prediction uncertainty.*Geophys. Res. Lett.***42**, 7949–7957 (2015). - 37.
Gusman, A. R.

*et al*. earthquake estimated from ocean-wide tsunami waveforms and GPS data: The 2014 Iquique earthquake.*Geophys. Res. Lett.***42**, 1053–1060 (2015). - 38.
Jaramillo, E.,

*et al*. Calibrating coseismic coastal land-level changes during the 2014 Iquique (Mw=8.2) earthquake (northern Chile) with leveling, GPS and intertidal biota.*PLoS One***12**, (2017). - 39.
Shrivastava, M. N.

*et al*. Coseismic slip and afterslip of the 2015 Mw 8.3 Illapel (Chile) earthquake determined from continuous GPS data.*Geophyiscal Res. Lett*. 710–719, https://doi.org/10.1002/2016GL070684 (2016). - 40.
Yue, H., Lay, T. & Koper, K. D. En échelon and orthogonal fault ruptures of the 11 April 2012 great intraplate earthquakes.

*Nature***490**, 245–249 (2012). - 41.
Satriano, C., Kiraly, E., Bernard, P. & Vilotte, J. The 2012 Mw 8. 6 Sumatra earthquake: Evidence of westward sequential seismic ruptures associated to the reactivation of a N-S ocean fabric.

*Geophyiscal Res. Lett*.**39**(2012). - 42.
Cahyadi, M. N. & Heki, K. Coseismic ionospheric disturbance of the large strike-slip earthquakes in North Sumatra in 2012: Mw dependence of the disturbance amplitudes.

*Geophys. J. Int.***200**, 116–129 (2015). - 43.
Wessel, P., Smith, W. H. F., Scharroo, R. & Luis, J. and F. W. Generic Mapping Tools: Improved Version Released.

*EOS Trans. AGU*94, (409–410 (2013). - 44.
Yagi, Y. & Okuwaki, R. Integrated seismic source model of the 2015 Gorkha, Nepal, earthquake.

*Geophys. Res. Lett.***42**, 6229–6235 (2015).

## Acknowledgements

A. S., D. T., and S. N. sincerely thank the Department of Science and Technology (DST), Govt. of India for the research fellowships. This work is part of the interdisciplinary program ‘Coupled Lithosphere-Atmosphere-Ionosphere-Magnetosphere System (CLAIMs)’ of Indian Institute of Geomagnetism and supported by DST, India. LR acknowledge the support of French national space agency (CNES) and Université Côte d’Azur (UCA).

## Author information

### Affiliations

### Contributions

M.S.B. conceived the idea of developing a 3D non-tectonic forcing mechanisms model along with L.R. M.S.B. designed and interpreted the results and drafted the text. L.R. contributed during the interpretation of results and provided useful suggestions. A.S. performed the analysis presented in the paper, prepared the figures and co-drafted the text. A.S., S.N., and M.P. carried out crustal deformation analysis. D.T. participated in scientific discussion. D.S.R. provided useful suggestions during the preparation of the results. All authors actively discussed the results and reviewed the manuscript.

### Corresponding author

## Ethics declarations

### Competing interests

The authors declare no competing interests.

## Additional information

**Publisher’s note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## Supporting Information

## Rights and permissions

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

## About this article

### Cite this article

Bagiya, M.S., Sunil, A.S., Rolland, L. *et al.* Mapping the Impact of Non-Tectonic Forcing mechanisms on GNSS measured Coseismic Ionospheric Perturbations.
*Sci Rep* **9, **18640 (2019). https://doi.org/10.1038/s41598-019-54354-0

Received:

Accepted:

Published:

## Further reading

## Comments

By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.