Total Electron Content (TEC) estimation over Pakistan from local GPS network using spherical harmonics

Many organizations allow GNSS users to access Global Ionosphere Maps (GIMS). However, the TEC estimates derived from GIMs are of insufficient quality to describe small-scale TEC variations over Pakistan. In this paper, the first local TEC map over Pakistan for the year 2019, derived from a regional GPS network, is presented. Spherical harmonics expansion is employed to estimate TEC with the spatial resolution of latitude 0.2° x longitude 0.2° and temporal resolution of 5 minutes. The impact of changing the degree/order of harmonics is assessed and it is determined that harmonic expansion up to 6 degrees is sufficient for estimating accurate TEC map for the region of interest. We have demonstrated that the TEC maps of Pakistan generated by local model conform better to the GIM by Center of Orbit Determination (CODE) ( RMS = 5.83 ) as compared to International Reference Ionosphere (IRI-2016) ( RMS = 7.18 ). We found that the TEC estimated by the local model shows a better correlation to measured TEC; CODE-GIM overestimated TEC, while IRI-2016 underestimates it. Moreover, it was observed that TEC peaks during noon (1100-0100 LT) and Equinox (April). The residuals of local TEC estimates with respect to TEC obtained from CODE-GIM indicate the inaccuracy of CODE-GIM over the region of Pakistan: highest deviation of TEC from local model with respect to CODE –GIM was observed in April ( RMS = 8.73 ) and minimum in October ( RMS = 2.78 ). We have also analyzed the performance of our maps in geomagnetically disturbed days. The research presented in this paper will contribute towards the ionosphere study over Pakistan, where limited research is available currently.


Introduction
Total Electron Content (TEC) from Global Navigation Satellite System (GNSS) has a number of significant applications in a variety of fields including but not limited to navigation, timing, surveying, agriculture, telecommunications, financial industry and military. As with all electromagnetic signals, GNSS signals are also prone to a number of natural and artificial error sources such as atmospheric errors, clock errors, multipath, hardware biases, jamming and spoofing. However, studies have confirmed the GNSS ionospheric delay as the single most influential cause of GNSS positioning error [Klobuchar, 1982;Feltens and Schaer, 1998;Krankowski et al. 2005;Perez, 2017;. Therefore, it becomes pertinent to model the ionosphere delay and variations and account for this error in positioning.
The ionosphere is a region of Earth's atmosphere containing charged particles along with free electrons having boundary limits starting from ~50 km and extending more than ~1000 km. This region causes a frequencydependent propagation effect on the RF signals passing through it, such as phase delay and carrier advance. One of the parameters to quantify the ionosphere is TEC which is the total number of free electrons present in a square meter cross sectional area between the satellite and the receiver.
The analysis of the received GNSS signal properties may contribute towards detection and identification of the ionospheric condition variations due to space weather events such as solar flares, solar radio bursts (SRBs) and coronal mass ejections (CME) [Cannon, 2013]. With the deployment of several GNSS stations around the globe, it has become an alternate source of data to model and monitor the ionosphere TEC. Moreover, the vast network of GNSS stations across the globe allows large scale spatial coverage thus resulting in a database enough to create global maps of TEC variability [Mendillo, 2006]. In fact, several international organizations use the GNSS data to provide global ionosphere maps (GIMs) representing the TEC distribution using below various techniques.
International GNSS Service (IGS) uses a weighted mean of other analysis centers to create a GIM in Bernese software . Jet propulsion laboratories (JPL) uses a three shell model based on spline functions. This approach uses the Kalman filter and 4DVar (four-dimensional variation) to estimate the state vectors and update the observation matrix [JPL, 2015]. CODE uses carrier-smoothed pseudo ranges in solar-geomagnetic reference frames and spherical harmonic expansion up to 15 degrees. European Space Agency (ESA) uses spherical harmonic functions expansion up to 15 degrees and assumes that the ionosphere is concentrated as a single-layer at 450km.
Technical University of Catalonia (UPC) uses a tomographic approach which involves feeding the GPS NTrip data streams from a global network of receivers to solve the main tomographic equation with the help of TOMION software: where S is the STEC, B is the unknown ambiguity, , is the mean electron density within each voxel and is the length intersection between each voxel.
These institutes estimate GIMs with a resolution of 5 o in latitude, 2.5 o in longitude and 2 hours in time [Hernández-Pajares et al., 2009]. However, the resolution of global TEC map is insufficient for derivation of the regional TEC map of required quality [Wielgosz et al., 2003]. On the other hand the empirical models, like IRI, are not efficient in describing certain ionosphere anomalies especially those that only occur locally.
A lot of effort has been made to enhance the knowledge of ionosphere variability at regional level by generation of local/regional TEC maps. The techniques for mapping TEC are mainly divided in to two categories: (1) Functionbased models and (2) Grid-based models. The former includes mathematical models such as spherical harmonics and B-splines while the latter approach divides the geographic region into different grids and interpolates data over each grid with no IPP. The regional TEC map over Canada has been developed using spherical harmonics from daily RINEX files recorded at an interval of 30 secs [Reza et al., 2011]. The regional TEC model over Africa referred to as SATCP (South African TEC prediction) has been developed using spherical harmonic technique which generates the model over Africa with a resolution of 1 minute and produces results comparable to that of ionosondes [Habarulema et al., 2011]. Neural Networks specifically multilayer perceptron (MLPs) and radial basis function networks (RBFN), have been used to demonstrate the TEC mapping over Turkey due to their non-linear modeling capability [Yilmaz et a., 2009]. The TEC variations over Japanese islands have been mapped with 0.5 o x 0.5 o grid using spherical harmonics expanded to as high as 60 degrees [Ping et al., 2002]. A recent study [Krypiak-Gregorczyk et al., 2017] used a novel approach utilizing thin splines to model the ionosphere TEC over Europe. The TEC over the region of China has been developed using tomography [Zhizhao and Gao, 2004] as well as spherical harmonics [Zhao et al., 2016]. Studies have also been carried out in Japan for estimating the TEC using spherical harmonics with GEONET (GNSS Earth Observation network) data [Ohashi et al., 2015].

Maria Mehmood et al.
The standard ionosphere models (e.g IRI and GIM) have poor performance in the equatorial or low latitude regions as well as regions with sparse data; which is why GNSS users in these regions cannot fully rely on standard models. This creates the need to have tools that will enable the local GNSS users with reliable TEC data. In this context, a local TEC model describing the variation of TEC for year 2019 is presented. The local TEC map has been generated with a spatial resolution of latitude 0.2° x longitude 0.2° and temporal resolution of 5 minutes, using spherical harmonic expansion up to 6 degrees. The focus of the research is not only to demonstrate the performance of standard models over Pakistan but to present an improved alternative in the form of locally-driven TEC maps.
These TEC maps provide an accurate assessment of the TEC variation over Pakistan for year 2019 as compared to the currently-used CODE-GIM. The main reason for this is because among all the IGS (International GNSS Service) stations used to generate the global map, not even one is located in Pakistan ( Figure 1). The results from the generated model are compared to CODE-GIM as well as IRI-2016 to determine the improvement offered in terms of TEC estimation.

Data set
The data set used for this research is collected from four GNSS stations located in Pakistan ( Figure 2). The selection of reference sites is such that the obtained Ionosphere Pierce Points (IPPs) cover the whole region of interest. Moreover, the monitoring equipment selection, installation and operation methods were in complaint with the guidelines provided by IGS [IGS, 2015]. The lack of IPPs over extreme south-west reveals that it will be advantageous to add another station towards the south-west of Pakistan in the future. The data is archived in RINEX format on daily basis with an interval of 30 seconds. The quality of logged data is ensured by passing it through TEQC software developed by UNAVCO, available at https://www.unavco.org/software/data-processing/teqc/teqc.html.

3
TEC estimation over Pakistan from local GPS using SH

TEC estimation from GPS
The RINEX file from a dual frequency GPS receiver includes carrier phase (Φ₁ and Φ₂) as well as pseudorange ( ₁ Where indicates the frequency, ₁ or ₂, ρ is the geometric range between the satellite and receiver in meters, is the orbital error in meters, are the satellite clock error and receiver clock error with respect to GPS time in seconds, is the troposphere error in meters, is the ionosphere delay for each frequency and are the satellite delay and receiver delays in meters for both carrier phase and pseudorange on each frequency, are the measurement noise errors in meters, are the integer ambiguities and is the speed of light (3.08 x 10 8 m/s).
The geometry-free linear combination of pseudo range is then obtained by estimating the difference between the pseudo ranges on both frequencies, as shown below: The geometry-free combination solution for TEC estimation yields the slant ionospheric TEC (STEC) which is expressed as: TEC is measured in terms of TECU where 1 TECU = 10 16 electrons/m2. The same procedure may be followed to obtain the TEC from carrier-phase measurements. In case of two-dimensional ionospheric modelling the STEC is converted to vertical TEC (VTEC) using the single layer model (SLM) (Figure 3).
SLM assumes that the ionosphere is concentrated as a single layer around 350-450 km above the Earth surface.
STEC is mapped on a vertical to obtain VTEC using an obliquity function at the ionosphere pierce point (IPP) defined by [Han et al., 2013]. The VTEC is expressed as: (6) where R is the equatorial radius of the earth (6378 km), H the altitude of the ionosphere shell and z the zenith angle [Davies, 1990]. For this study the STEC was mapped onto VTEC using SLM with a height of 350 km.
The spherical harmonic expansion is applied to VTEC as described in [Zhao et al., 2016]: where and are the unknown model co-efficients. θ and ξ are the magnetic latitude and local time ionosphere pierce points respectively in the spherical cap harmonic coordinate system, which are obtained from the measured longitude and latitude of ionosphere pierce points. The degree and order of spherical harmonics are represented by n and m and signify the resolution in latitude and longitude respectively. is the normalized associated Legendre function and maybe expressed as = Λ (n, m) where Λ refers to the normalization function being applied on the un-normalized polynomial . The normalization function may be expanded as: where is the Kronecker Delta. Substituting the equation 5 & 6 into 7, results in: The four unknown parameters: , , and , are estimated using least square (LS) method. An external constraint is assumed that the sum of all GPS satellite DCBs is zero. This is known as the zero-mean constraint and separates the receiver DCBs from satellite DCBs as well as makes Eq. 9 full rank.

TEC estimation from CODE
The CODE of IGS analysis center is operated by the Astronomical Institute of the University of Bern (AIUB), Switzerland. CODE uses data from more than 200 GNSS stations ( Figure 1) and uploads daily Global Ionosphere Maps (GIM) on their FTP http://ftp.aiub.ch/CODE/IONO/. The ionosphere data is disseminated in the form of TEC values in IONEX [Hernández-Pajares et al., 2009] or GTEX [Saito et al., 2017] format. A detailed description of these products is available on their website at http://cmslive3.unibe.ch/unibe/philnat/aiub/content/e15/ e59/e440/e447/e458/index_eng.html. The TEC data in the IONEX files is sorted according to the geographical latitude and longitude. The latitudes range from -180° to 180° with a step size of 5° and the longitude range from -87.5° to 87.5° with a step size of 2.5°; this makes a total of 5183 grid points (each with a unique estimated TEC value) in a CODE GIM. These values are then interpolated to estimate TEC values all over the globe.
However, since there are no grid points (or IGS stations) with in Pakistan (Figure 1) any TEC estimation over this region using CODE-GIM will have significant interpolation error. Poor spatial coverage and large deviations from a regional model, where less or no local GNSS data was used were also observed by [Ohashi et al., 2015]. For this study the IONEX files for the year 2019 were download from the aforementioned website and interpolated over the region of interest via post-processing.
The IRI model divides the ionosphere into six sub-regions and provides information on several ionosphere parameters including electron density and temperature, ion composition, drift and temperature and TEC up to a height of 2000 km. The IRI model is the most-widely used empirical and deterministic model to observe ionosphere trend for space weather. The advantage of this model is that since it is empirical, it is independent of our changing knowledge about the physics of atmosphere. The disadvantage lies in the fact that its accuracy depends on its database, so the geographical and temporal locations not covered by the database will experience error. The database of IRI is concentrated in mid latitudes, therefore low-latitude regions experience error [Bilitza, 2018]. The latest IRI model is IRI-2016 and a detailed description is available at http://irimodel.org/.

Results and discussion
This study presents a local TEC map for Pakistan (20° -40° latitude, 55° -80° latitude), as shown in Figure 2, for year 2019, using spherical harmonic functions. The local map provides a resolution of 0.2° x 0.2° x 5 mins and is estimated at a height of 350 km. For global representation of TEC the degree and order are set higher than 10 for generation of global map with sufficient resolution, but for regional mapping they need to be determined accordingly [Zhao et al., 2016]. For this study the spherical harmonic expansion of several degrees was examined starting from 3 and going upwards. The expansion of up to 6 satisfies the previously set requirements for the experimental local TEC map accuracy over Pakistan, since there is no significant change in RMS error with respect to CODE-GIM after further expansion (Table 1).
Maria Mehmood et al.  GIM as compared to those by IRI-2016; Moreover, the difference between this study and CODE-GIM is lesser than that between CODE-GIM and IRI-2016. Similar results were observed by [Xin et al., 2016]   The same strategy was extended to obtain the residuals between CODE-This study and IRI-This study for every hour of each day of the year 2019. The histogram of the residuals ( Figure 5), with extreme outliers ignored, reveals that the residuals between CODE and this study are positive with a mean error of 7.26 and variance of 3.31; whereas the residuals of IRI and this study were negative with a mean of -11.99 and variance of 2.41. The residual analysis confirms that the results obtained from this study conform better to CODE GIM rather than IRI-2016 and that CODE GIM overestimates the TEC whereas IRI-2016 underestimates it.
In order to validate the model with actual data, a GPS scintillation receiver was set up at a test location of latitude   (Figure 7). The TEC residuals at each hour is calculated by taking the difference between the two hourly means of both models at same hour thus: where = hour number.
It can be observed from Figure 7, that the residuals are higher during the day as compared to night and maximum during Equinoxes. Furthermore, the daily mean of TEC estimated by this study and CODE-GIM at the test location was observed and their difference was calculated (Figure 8). The residual TEC was then divided into months and a monthly RMS of the residual TEC was calculated. A maximum RMS of 8.736 TECU was recorded in the Equinox (April) and minimum of 2.780 TECU was obtained during winter solstice (October).

9
TEC estimation over Pakistan from local GPS using SH  The VTEC maps were generated using this study over 24 hours for DOY 312 (Year 2019) to assess the diurnal variation ( Figure 9). It can be see that VTEC values increase slowly in the morning with sunrise and decrease with sunset in the evening. The maxima is obtained around noon (1100 LT-0100 LT), which may be due to the intensification of photoionization process as described in [Chapman, 1931]. A similar pattern was observed by previous studies over the region where VTEC was estimated using GPS and IRI (International Reference Ionosphere) [Tariq et al., 2019].
In order to emphasize the importance of having a local TEC model for our region under consideration, we tested

Conclusion
Currently, TEC maps provided by IRI-2016 and CODE-GIM are being used for TEC variations and estimation within Pakistan. However, these models not only offer a low resolution in space and time but exhibit poor performance since no GNSS station is located within Pakistan. Due to lack of GNSS-related research and infrastructure in the region of interest no customized solution for TEC maps is currently available for this area.
This study has developed a local TEC map with a resolution of latitude 0.2° x longitude 0.2° x 5 mins, for the first time, over the region of Pakistan. The research has also determined that for GPS-based TEC mapping of this region the spherical harmonic function may be expanded up to 6 degrees for optimal results. The TEC values estimated by the local TEC maps are compared to the presently used TEC maps (provided by CODE-GIM and IRI-2016) as well as actual TEC measurements from a GPS scintillation receiver. It was found that the local TEC maps by this study are an accurate estimation of TEC for this region for Year 2019 as compared to CODE-GIM and IRI-2016. The study also found that the local TEC maps correlate better to CODE-GIM rather than IRI-2016. We have also analyzed the performance of all three models in case of a G3 geostorm that occurred on 14 th May 2019 and concluded that the local model exhibits much pronounced enhancements in TEC as compared to standard models.
The global models are focused on the universal distribution and climatological analysis of TEC and result in missing out small-scale local variations due to sparse distribution of stations. We have demonstrated that the regional maps are able to reproduce local variability, especially around the EIA. The TEC maps presented in this study may be utilized to accurately estimate the ionosphere delay for GPS users in Pakistan and improve the positioning solutions as well as to observe the geo storms that occur over this region. The study will contribute towards ionosphere mapping and forecasting service over Pakistan by allowing reproduction of local ionosphere morphology.
For future work, we recommend generation of local TEC maps using multi-constellation data (e.g. Galileo/ BDS). It is also suggested to explore the advantage of TEC calculation with geostationary satellite data [Savastano et al., 2019].