The natural hazard such as earthquakes, tsunamis and volcanic eruptions can cause severe loss of life and property each year. Therefore, rapid and accurate detection of natural hazard is very important to reduce its effect. To achieve this goal, multiple instruments and techniques are necessary. In recent year, a new natural hazard detection method, namely ionospheric sounding, attracts extensive attention of scientific communities. The natural hazard-related atmospheric acoustic waves or gravity waves can propagate to the ionosphere and interact with the plasma, then generate the traveling ionospheric disturbances (TIDs) (Hines 1972). When the electromagnetic waves such as Global Navigation Satellite System (GNSS) signals travel through the ionosphere, the TIDs can affect the electromagnetic waves and produce measureable values with it.

The first natural hazard-induced TIDs were reported after the Alaskan earthquake of March 28, 1964 observed by Doppler sounders and ionsondes (Davies and Baker 1965). In last two decades, owing to the high spatiotemporal resolution, GNSS total electron content (TEC) has widely been used to detect various natural hazards, including the earthquakes (Heki and Ping 2005; Astafyeva et al. 2009; Liu et al. 2011; Sanchez et al. 2022), tsunamis (Artru et al. 2005; DasGupta et al. 2006; Rolland et al. 2010; Tang et al. 2018) and volcanic eruptions (Heki 2006; Shults et al. 2016; Cahyadi et al. 2020; Themens et al. 2022). These observed results show the natural hazard signatures in ionosphere have similar propagation characteristics in term of period, horizontal speed and orientation compared to the natural hazard-induced acoustic waves or gravity waves. In addition, simulated results on natural hazard-ionosphere coupling process can further provide convincing evidence that the ionosphere is sensitive to the natural hazards (Occhipinti et al. 2008; Rolland et al. 2011).

The observed and simulated results show ionosphere sounding by GNSS TEC is a powerful tool to detect the natural hazard, which can serve as an important decision support for hazard warning and damage mitigation. Currently, there are many open-source programs for GNSS positioning and related applications (Nievinski and Larson 2014; Bahadur and Nohutcu 2018; Zhao et al. 2021; Altuntas and Tunalioglu 2022). However, open-source program on detection of natural hazard by GNSS TEC is very scarce. In view of this, the IonKit-NH toolkit is developed here.

### Ionosphere Sounding By Gnss Tec

Ionospheric TEC is the integrated electron density along the ray path from a satellite to a receiver. The slant ionospheric TEC can be computed using the dual-frequency GNSS observations for each satellite-receiver pair:

$$TEC=\frac{1}{A}\frac{{f}_{1}^{2}{f}_{2}^{2}}{{f}_{1}^{2}-{f}_{2}^{2}}\left({L}_{1}-{L}_{2}+N+\epsilon \right)$$

1

where \(A=40.3\) m3/s2, \({f}_{1}\) and \({f}_{2}\) are the carrier phase frequencies, \({L}_{1}\) and \({L}_{2}\) are carrier phase observations in meters, \(N\) is the unknown carrier phase ambiguity in meters, and \(\epsilon\) is the measurement noise. The carrier phase ambiguity is constant and eliminable in the following detrend process. To enhance usability, the approximate value of \(N\) is obtained using following equation

$$N=⟨{L}_{1}-{L}_{2}-\left({P}_{2}-{P}_{1}\right)⟩$$

2

where the notation \(⟨\bullet ⟩\) means the average operation, \({P}_{1}\) and \({P}_{2}\) are code observations.

The unit of TEC is TECU, which is equal to 1016/m2. The ionosphere is usually assumed as a thin shell layer (see Fig. 1). Then, the so called ionospheric pierce point (IPP), which is the intersection point between the ray path from a satellite to a receiver and the ionosphere shell, can be applied to express the position of TEC. The coordinate of IPP can be computed using following equation

$$\left\{\begin{array}{c}{\phi }_{p}={\text{sin}}^{-1}(\text{sin}{\phi }_{0}\text{cos}\varDelta z+\text{cos}{\phi }_{0}\text{sin}\varDelta z\text{cos}{A}_{z})\\ {\lambda }_{p}={\lambda }_{0}+{\text{sin}}^{-1}\left(\frac{\text{sin}{\Delta }z\text{sin}{A}_{z}}{{cos}{\phi }_{p}}\right)\\ \varDelta z=z-{{sin}}^{-1}\left(\frac{R}{R+H}{sin}z\right)\end{array}\right.$$

3

Where \({\phi }_{p}\) and \({\lambda }_{p}\) are the latitude and longitude of IPP, \({\phi }_{0}\) and \({\lambda }_{0}\) are the latitude and longitude of receiver, \(z\) and \({A}_{z}\) are the zenith angle and azimuth angle of the ray path at receiver, \(R\) is the earth radius and \(H\) is the height of ionosphere shell.

To extract the TIDs signal in ionosphere, the TEC variation (dTEC) time series is computed by the following equation

where \({TEC}_{f}\left(t\right)\) is the fitted TEC by the Savitzky-Golay smoothing filter (Savitzky and Golay 1964). After obtaining the dTEC time series from all the satellite-receiver pairs, the two-dimensional map and time-distance map of dTEC could be plotted. The dTEC two-dimensional map is plotted using the all the IPPs at a specific epoch, which is convenient to observe the propagation direction and distribution of the TIDs. As for the time-distance map, the dTEC is plotted as a function of distance from event center and time. Then, the slope of the line along the propagation direction of TIDs in time-distance map is the horizontal velocity.

### Ionkit-nh Toolkit

IonKit-NH toolkit is developed in MATLAB R2010b version, which has a wide range of compatibility. The main flowchart of ionospheric detection of natural hazard using IonKit-NH toolkit is shown in Fig. 2. Specifically,

*Step 1: Input GNSS files and set processing options*. The input files include observation and navigation files. Both RINEX version 2 and RINEX version 3 files are supported. The processing options include the type of dual-frequency combination, the height of ionosphere shell, the navigation system, the elevation mask and the thresholds of cycle-slip detection.

*Step 2: Read GNSS data*. The multi-system GNSS including GPS, GLONASS, Galileo and BDS are supported. The main functions include *readobs*, *setobstype* (called by *readobs*) and *readnav*.

*Step 3: Generate TEC files*. This is achieved by the function *calTEC*. In addition, three functions, namely *getfreqwave*, *getippinf* and *preprocess* are called. The function *getfreqwave* is to acquire the frequency and wavelength for each carrier wave signal. The function *getippinf* is to compute satellite elevation and azimuth and information on IPP including project factor, latitude and longitude. To get the satellite position, it further calls the functions *satpos* (for GPS, Galileo and BDS) or *satposg* (for GLONASS). The function *preprocess* is to check data validity and detect cycle slips.

*Step 4: Generate dTEC files*. This is achieved by the function *sgdTEC*. The event position (longitude and latitude) is necessary to calculate the horizontal distances between event center and IPPs. The order and sliding window of Savitzky-Golay smoothing filter can be set according to the event type.

*Step 5: Plot two-dimensional maps and time-distance maps of dTEC*. These are finished by the functions *ddmap* and *tdmap*, respectively. Here, the M_Map toolbox developed by Pawlowicz is called by function *ddmap*. The detailed instruction for the M_Map toolbox is available online at www.eoas.ubc.ca/~rich/map.html.

*Step 6: Obtain the propagation characteristics of TIDs*. Using the dTEC two-dimensional maps and time-distance maps, the propagation characteristics of TIDs are obtained. According to the propagation characteristics, we can determine if the TIDs are related to the natural hazard event.

For convenience, the script *getTEC* that integrates the Step 1–3 is provided in IonKit-NH toolkit. In addition, it also provides a script *plot2dmap* to make an animation (including GIF and video) of the dTEC two-dimensional maps, which is more vivid to show the TIDs’ evolutions. The users can change the options in the two scripts to obtain the results their need.