# Global GNSS-RO Electron Density in the Lower Ionosphere
## 2 Data and Methods
### \(N_{e}\) Retrievals  
#### 2.3.1 Inversion from \(\Delta\)hTEC  
The inversion algorithm for the \(N_{e}\) retrieval in this study uses the optimal estimation inversion method [45]. The measurement vector \(\mathbf{y}\) is the \(\Delta hTEC\) derived from Equations (3)-(8) as a function of \(h_{t}\), i.e., \(\mathbf{y}\equiv\{\Delta hTEC(h_{t});\,h_{t}=h_{1}\cdots h_{max}\}\), where \(h_{max}\) is the maximum \(h_{t}\) in the RO \(\phi_{ex}\) profile. The state vector \(\mathbf{x}\) is \(N_{e}\) profile, i.e., \(\mathbf{x}\equiv\{N_{e}(z);\,z=60\,\,\text{km},\,\cdots,z_{max}\}\). When the measurement vector is changed from \(hTEC\) to \(\Delta hTEC\), the weighting functions of the bottom-up algorithm are also transformed from the Abel functions to a new set of functions, as shown in red in Figure 3d. Unlike the Abel function, the new weighting functions for \(\Delta hTEC\) have a very small contribution from higher altitudes, which provide the advantage of minimizing \(F\)-region impacts. The linearized forward model for the measurement vector (\(\mathbf{y}\)) is related to the state vector (\(\mathbf{x}\)) through a weighting function \(\mathbf{K}\), or the _Jacobians_ as in \(\mathbf{K}=\partial\mathbf{y}/\partial\mathbf{x}\), as follows:  
\[\mathbf{y}=\mathbf{y}_{0}+\mathbf{K}\mathbf{\cdot}(\mathbf{x}-\mathbf{x}_{0} )+\epsilon_{\mathbf{y}} \tag{10}\]  
Figure 3: (**a**) State vector \(\mathbf{x}_{0}\), (**b**) measurement vector (\(\Delta hTEC=hTEC-hTEC_{fit}\) ), (**c**) \(\mathbf{y}_{0}\) in the linearization model, and (**d**) the associated weighting functions \(\mathbf{K}\) (red) and their comparison with the Abel functions (black) for selected tangent heights between 50 and 500 km. The weighting functions \(\mathbf{K}\) are flipped in sign purposely, to help the comparison. The intermediate step for the calculated RO \(hTEC\) and its linear fit (red line) from \(h_{t}=30\)–60 km is illustrated in (**b**), from which \(\mathbf{y}_{0}\) is derived as their difference. See text for more discussions. where the linearization is based on the state vector \(\mathbf{x}_{0}\) (\(N_{e}\) in m\({}^{-3}\)) and the measurement vector \(\mathbf{y}_{0}\) (\(\Delta hTEC\) in TECu). A single \(N_{e}\) profile is used for the \(\mathbf{x}_{0}\) in Equation (10), which is a merged profile from the annual mean of IRI-2016 from 2008 for the \(F\)-region [30] and the mean FIRI \(E\)-region profile for the solar minimum condition [29]. We assume that the weighting function \(\mathbf{K}\) is invariant with the coefficient \(a_{1}\) (i.e., the slope of the linear function in Equation (8)). The linearization model \(\mathbf{x}_{0}\), \(\mathbf{y}_{0}\), and \(\mathbf{K}\) are depicted in Figure 3. The measurement uncertainty of \(\mathbf{y}\), or \(\varepsilon_{y}\), is 0.3 m for the \(\phi_{exL1}\), or 1.9 TECu for \(\Delta hTEC\), to account for multipath, RIE, scintillation, neglection of ionospheric bending, and other measurement errors. The multipath error from reflections off the spacecraft structure is estimated to <0.2 m for COSMIC-2 [46] and <1 m for COSMIC-1 [38]. The RIE error, as shown in Figure 1, is on the order of \(\sim\)0.1 m but can be larger in a strong scintillation condition. The \(\mathbf{x}\) inversion in the optimal estimation method [45] can be expressed by  
\[\hat{\mathbf{x}}=\left[\mathbf{S}_{\mathbf{a}}^{-1}+\mathbf{K}^{T}\mathbf{S}_ {y}^{-1}\mathbf{K}\right]^{-1}\left[\mathbf{S}_{\mathbf{a}}^{-1}\mathbf{a}+ \mathbf{K}^{T}\mathbf{S}_{y}^{-1}\mathbf{y}\right] \tag{11}\]  
where \(\hat{\mathbf{x}}\) is an optimal solution to \(\mathbf{x}\) in Equation 10 and \(\mathbf{a}\) is the a priori of \(\mathbf{x}\), which is set to be same as \(\mathbf{x}_{0}\) in this study. \(\mathbf{S}_{\mathbf{a}}=0.02\cdot(5\cdot 10^{8}+\mathbf{x}_{0})\cdot I\) and \(\mathbf{S}_{y}=\varepsilon_{y}^{2}\cdot I\) are covariance matrices for \(\mathbf{a}\) and \(\mathbf{y}\), respectively.