# Monitoring of Ionospheric Anomalies Using GNSS Observations to Detect Earthquake Precursors
## 2 Materials and Methods
### Model Implementation: From GNSS Doppler Observations to VTEC and VTEC Rate
#### 2.1.1 Equations of GNSS Doppler Observations  
A GNSS Doppler measurement is the rate of the phase observation. So, if we use \(j\) to identify carrier frequencies (\(j=1,2\)), \(r\) receivers, and \(s\) satellites, the Doppler observation at the time \(t\) is \(D^{s}_{r,j}(t)=-\lambda_{j}\Phi_{j}(t)\)[3; 46; 47], where \(D^{s}_{r,j}(t)\) is the Doppler observation, \(\Phi_{j}(t)\) is the phase time derivative, and \(\lambda_{j}\) is the carrier wavelength. The corresponding observation equation is as follows:  
\[D^{s}_{r,j}(t)=\dot{\rho}^{s}_{r,j}(t)+c\dot{\rho}^{s}_{r}(t)-\dot{l}^{s}_{r,j} (t)+\dot{l}^{s}_{r,j}(t)+c\Big{(}\dot{\varepsilon}^{s}_{j}(t)+\dot{\varepsilon }_{r,j}(t)\Big{)}+\dot{\varepsilon}^{s}_{r,j}(t). \tag{1}\]  
In Equation (1), \(\dot{\rho}^{s}_{r,j}(t)\) is the velocity of the range between the satellite and receiver, \(\dot{\delta}^{s}_{r}(t)\) is the rate of the satellite-receiver clock error combination, \(\dot{l}^{s}_{r,j}(t)\) and \(\dot{l}^{s}_{r,j}(t)\) are the rates of the ionospheric and tropospheric terms, \(\dot{\varepsilon}^{s}_{j}(t)\) and \(\dot{\varepsilon}_{r,j}(t)\) are the rates of frequency-dependent satellite and receiver biases [48], and \(\dot{\varepsilon}^{s}_{r,j}(t)\) is the noise term; all the terms are in m/s. The ionospheric error can be written as a function of VTEC and the zenith angle at the Ionosphere Pierce Point (IPP) \(z^{\prime}:\;\dot{l}^{s}_{r,j}(t)=\frac{1}{\cos z^{\prime s}_{r}(t)}\frac{40.3 }{f_{j}^{2}}VTEC(t)\), where \(f_{j}\) is the broadcasting carrier frequency in Hz, and VTEC is the vertical TEC in TECU (1 TECU = 10\({}^{16}\) electrons per m\({}^{2}\).). The zenith angle depends on the height of the ionospheric layer \(h_{m}\), the Earth radius \(R_{e}\), and the zenith angle at the observation site \(z^{s}_{0r}(t)\): \(\sin z^{\prime s}_{r}(t)=\frac{R_{e}}{R_{e}+h_{m}}\sin z^{s}_{0r}(t)\)[46]. The corresponding time derivative is as follows:  
\[\dot{l}^{s}_{r,j}(t)=\frac{\dot{z}^{\prime s}_{r}(t)\tan\big{(}z^{\prime s}_ {r}(t)\big{)}}{\cos z^{\prime s}_{r}(t)}\frac{40.3}{f_{j}^{2}}VTEC(t)+\frac{1 }{\cos z^{\prime s}_{r}(t)}\frac{40.3}{f_{j}^{2}}VTEC(t), \tag{2}\]  
where \(\dot{z}^{\prime s}_{r}(t)=\frac{R_{e}}{R_{e}+h_{m}}\frac{\cos z^{s}_{0r}(t)} {\cos z^{\prime s}_{r}(t)}\dot{\varepsilon}^{s}_{0r}(t)\) and \(\dot{z}^{s}_{0r}(t)=\frac{\dot{\varepsilon}^{s}_{0r}(t)\cos z^{s}_{0r}(t)- \dot{u}^{s}_{r}(t)}{r^{s}_{r}(t)\sin z^{s}_{up component of \(r_{r}^{s}(t)\) is then \(u_{r}^{s}(t)=r_{r}^{s}(t)\cos z_{0r}^{s}(t)\), where \(r_{r}^{s}(t)\), \(r_{r}^{s}(t)\), \(u_{r}^{s}(t)\), and \(u_{r}^{s}(t)\) are functions of the satellite position in the Earth-Centered, Earth-Fixed (ECEF) reference system. If we use broadcast orbits, the rates can be computed as described in [PERSON] (2016); if instead we start from precise orbits, Lagrange interpolation can be used to compute precise coordinates and then the corresponding rates via a numerical differentiation method (e.g., using a five-point formula ([PERSON], 2017)). To generate the geometry-free combination, we subtract L1 and L2 Doppler observations \(D_{r,A}^{s}(t)=D_{r,1}^{s}(t)-D_{r,2}^{s}(t)\), and we introduce the rates of differential code biases ([PERSON], 2017; [PERSON], 2017) (inter-frequency biases) for the satellite and for the receiver, respectively, \(D\hat{C}B^{s}(t)=\varepsilon_{1}^{s}(t)-\varepsilon_{2}^{s}(t)\) and \(D\hat{C}B_{r}(t)=\varepsilon_{r,1}(t)-\varepsilon_{r,2}(t)\). At this stage, the observation equation is finally as follows:  
\[\begin{split} D_{r,A}^{s}(t)=&\frac{40.3(f_{1}^{2}- f_{2}^{2})}{f_{1}^{2}f_{2}^{2}\cos z_{r}^{s,s}(t)}z_{r}^{s}(t)\tan z_{r}^{s}(t) VTEC(t)+\\ &+\frac{40.3(f_{1}^{2}-f_{2}^{2})}{f_{1}^{2}f_{2}^{2}\cos z_{r}^ {s,s}(t)}VTEC(t)+c\left(D\hat{C}B_{r}(t)-D\hat{C}B^{s}(t)\right)+e_{r,A}^{s}(t).\end{split} \tag{3}\]