A Novel Approach to Study Regional Ionospheric Variations Using a Real-Time TEC Model

Since IGY (International Geophysical Year), through coordinated global observations, ionospheric research has been carried out by many countries. This effort primarily helped in the design and operation of HF radio wave communication systems. The Indian region covers a highly variable part of the equatorial electrojet and EIA (Equatorial Ionisation Anomaly) phenomena making its predictability difficult. With the advent of satellite communication and navigation, the need for accurate ionospheric TEC (Total Electron Content) models at global and regional scales has been stressed. The GAGAN (GPS Aided Geo Augmented Navigation) project jointly undertaken by the Indian Space Research Organisation (ISRO) and the Airport Authority of India (AAI) aims at effectively utilising the Global Navigational Satellite System (GNSS) to determine position coordinates accurately for aircraft precision landing applications. For this purpose the range errors are estimated by using a ground network of TEC stations spread over Indian region. The near simultaneous data collected from these dual frequency GPS stations can be used to generate the geo-referenced TEC values for various applications. The author has developed necessary algorithm and associated computer programmes for a real-time vertical TEC (VTEC) model based on TEC data collected from the GAGAN ground based network stations. The model has been tested and sample results presented here show that it adequately provides for the latitudinal resolution of 1 ̊ for the entire longitude span and also for two longitude blocks (73 83 & 83 93 ̊E) separately. Cubic spline and bilinear interpolation techniques are used for filling up temporal and spatial data gaps. The model provides tabulated output of hourly average VTEC data with latitude for ready use, as well as graphical displays of VTEC maps and contours for monitoring purpose. The real-time model and its extensions are also being used for detailed scientific studies; examples of these show small day to day variability of VTEC without any change in solar activity and indication of the change in the shape of the VTEC diurnal curve with season. The present model will be used for further studies to derive the monthly average variation of the diurnal pattern and the relationship between VTEC peak amplitudes with changes in solar activity. The new information generated can be fed back to improve the real-time model so that eventually the dependence of such models on ground based network stations data can be minimised.


Introduction
The region of earth's upper atmosphere covering ~50 -1000 km where free electrons exist was named "ionosphere" by Watson Watt in 1926 [1].The ionosphere is formed mainly by the interaction of solar extreme ultra violet (EUV) and X-ray radiations with the neutral atmospheric constituents like atomic oxygen (O), nitric oxide (NO), molecular oxygen (O 2 ) etc.The topmost F 2 layer of the ionosphere contains the bulk of the height integrated electron density and reaches the daytime peak values (~10 6 cm −3 ) around 350 km altitude [2].During night time, the F 2 layer also persists with reduction in the electron densities due to (a) the absence of solar ionising radiation and (b) continuation of night time ion-chemical loss processes at a slower rate [3].Prior to the advent of satellite based communication technology, the property of the F 2 layer to reflect radio waves was effectively utilised for 24-hour HF (3 -30 MHz) radio broadcasting and communication applications all over the world.Except under some special conditions of ducted propagation enabled through the presence of E and F region electron density irregularities caused by plasma turbulence [4], radio waves having frequencies greater than ~25 -30 MHz penetrate through the ionosphere and this property is utilised for trans-ionospheric satellite communication.
Unlike the visible and infrared radiations received on the earth, the solar ionising radiations (absorbed in the upper atmosphere and by the stratospheric ozone layer) are highly variable depending on both long term (11-year solar cycle) and short term energetic phenomena (solar flares, coronal mass ejections, etc.) taking place on the sun [5].The solar cycle variations change the electronion production rates thus impacting on the electron density profiles of the ionosphere by a factor of 5 -6 and even more if there are energetic short period events [2].The ionospheric plasma also varies due to the electrodynamical effects of ambient and imposed electric and magnetic fields.The scientific investigations of ionosphere have been mainly aided by developing ground based radio probing and radio propagation techniques and sound-ing rocket in-situ measurement devices to understand its layered regular and irregular structures, diurnal, seasonal, latitudinal and solar activity related variations.Due to the special horizontal configuration of the earth's magnetic field near the geographic equator, the low latitude ionosphere shows large variability due to the equatorial electrojet [6] and the EIA [7,8] related phenomena.
For ensuring the reliability of HF radio broadcasting and communication, it was necessary to undertake extensive ionospheric measurements at global level.Due to the more variable nature of the ionosphere over low and polar latitudes, a major emphasis was given to improve our understanding of its morphology and various physical phenomena in a coordinated manner.In view of this important requirement, a world-wide IGY programme was launched as the first large scale global coordinated effort of observations using similar instruments during 1957-58.Networks of ground based radio probing techniques like ionosonde, ionospheric absorption instrument, riometer, airglow technique etc., were employed, each participating country taking the responsibility to operate, analyse, study and deposit data to the World Data Centre (WDC) [9]. Figure 1 shows the Indian activities related to equatorial and low latitude ionospheric research since IGY superimposed on the plot for the solar 10.7 cm radio flux over the long stretch from IGY to the recent times showing a number of 11-year solar cycle activity periods.The progress from ground based radio probing techniques to use of rocket and satellite platforms with participation in coordinated national and international programmes have been the main tenet of the Indian activity in this field.The shifting of emphasis from HF radio to microwave satellite broadcasting/communication temporarily reduced the focus on ionospheric research but soon it was realised that even for the trans-ionospheric microwave signals, the assessment of the range errors due to its passage through the ionosphere needs more accurate knowledge of the medium particularly for navigational satellite applications.The problem gets compounded not only due to the temporal and spatial variations of TEC causing unpredictable signal phase delays but also the presence of ionisation irregularities producing fading/loss of signal thus impairing the satellite navigational technique.

Satellite Navigational Systems
The GNSS being employed by many countries for regional and global applications (like the GPS) have their inbuilt system to apply for range corrections (i.e. for single frequency users) through a network of TEC/reference stations.The dual frequency (L1:1575.42MHz & L2: 1227.60 MHz) reception of GPS helps to derive accurate TEC values (1 TEC unit = 10 16 electrons/m 2 ) which can be used to generate a real time model of TEC [17] in addition to using a number of available TEC models from past and continuing efforts [18].This homogenised observation system has also served as a unique opportunity to study the high resolution ionospheric phenomena using the multi-station TEC data sets, the type of which have only become available during the GPS era.
In any Satellite Based Augmentation System (SBAS), the navigational satellite signal range errors are estimated from the line-of-sight or slant ionospheric delay measurements (1 TEC unit equals 0.16 m in range error at L1 frequency of GPS) by ground network stations in real time and broadcast to users, e.g.pilots for accurate position finding en route and for precision landing.The dual frequency pseudo range and carrier phase measurements for each station are converted into range errors pertaining to a sequence of ionospheric pierce points covered by a navigational satellite pass.The ionospheric correction information is provided at the corners of the matrix of 5˚ × 5˚ latitude-longitude grid points.The user can estimate his 3D coordinates using a single frequency receiver data and by converting the Vertical TEC (VTEC) values broadcast by SBAS (through a geostationary satellite) to Slant TEC (STEC) values using a suitable mapping function [19].
In addition to presenting the results of the model grid point values of VTEC of some selected days, the present work also demonstrates a new approach to investigate the regional TEC variations with high time-space resolutions using extended algorithm of the real-time model already developed by the author [20] for the Indian low latitude region and measured TEC data for the period May, 2007-April, 2008 from the dual frequency GPS receiver network stations obtained through the Indian GAGAN project.

GAGAN Project, Real-Time Algorithm/Model, Data and Method of Analysis
GAGAN is a joint project of ISRO and the AAI to demonstrate the Wide Area Augmentation System (WAAS) over Indian airspace.The main aim of this project is to provide Category-I Positioning Accuracy and Integrity to all the civilian aircrafts, in Indian Airspace, while landing.The requirement of accuracy in correcting the single frequency GPS receiver range error called Grid Ionosphere Vertical Delay (GIVD) should be <0.5 m (or TEC < 4 TEC units) for such applications.A user/pilot of an aircraft, depending on the grid in which he is located, uses the GIVD at the corners of the grid to determine the ionospheric line-of-sight delay with reference to the coordinates of the Ionospheric Pierce Point (IPP).Due to sharp and highly variable latitudinal gradients of TEC over the Indian low latitude region it is not possible to obtain accurate TEC values from the existing low resolution global models.Hence GAGAN approach is to develop a real time TEC model based on continuous observations of dual frequency GPS signals over a network of reference ground stations [16].Each network receiver collects raw pseudo range and carrier phase measurements round the clock from all the visible GPS satellites above >15˚ elevation angle.The raw measurements are converted to STEC at intervals of one minute and stored in a PC.This data is archived in a server for further data analysis.In the Technology Development phase of GA-GAN, 18 stations have been operating for TEC and scintillation measurements since April, 2004.The details of data processing and archival in a standard format are being carried out by the project [21].The effort to develop a 5˚ × 5˚ latitude-longitude grid based real-time TEC model by GAGAN included making some of the archived data of the individual TEC network stations available to scientists/academicians for detailed investigations.Using such an opportunity of the available archived GAGAN data for the period May 2007-April 2008, the author developed the necessary algorithm for a real-time model and.modified it for carrying out studies with this unique data set.The map in Figure 2 shows the distribution of TEC ground stations over Indian region under the GAGAN project.
From one location, a TEC receiver can observe at least 8 GPS satellites and measure STECs in 8 different directions.These STEC values can be converted to VTEC values which correspond to the latitude and longitude at the IPPs which will change as the satellite moves in its orbit.The main ground network instrument consist of dual frequency pseudo range and carrier phase measuring P-code GPS receivers recording TEC and ionospheric scintillation indices in standard format.
The dual frequency GPS receivers can record either raw data every one second from all the simultaneously visible GPS satellites and/or a compressed data (e.g.every minute).The raw data collection rate is approximately 400 bytes, per satellite, per second.However, for TEC modelling and ionospheric scintillation studies, measurements taken every minute are considered adequate and hence the raw data can be recorded every minute.The data recorded include parameters like, week no., time of week in seconds, satellite no.(PRN), satellite azimuth and elevation, S 4 scintillation index, STEC values at 45, 30, 15 and 60 sec (now) in TEC units with differentials of TEC every 15 seconds, and signal lock times.The total data volume is about 2 M Bytes per day per station.These daily files per station are archived and transmitted to the master control station where a grid based ionospheric model can be run using available data from network stations.Uncorrected positions determined from GPS satellite signals produce range errors between 0 -50 meters and the aim is to reduce this error using the realtime TEC model results up to ~0 -0.5 metres.
In the context of generating a real-time TEC model, the author has developed the necessary algorithm and MATLAB codes which provide options for selecting satellite elevation angles above a certain value (the archiv-ed data includes all satellite passes above 15˚ satellite elevation angles) and desired time averaging (default as hourly averages at hourly interval).The time data gaps are filled up using cubic spline fitting for individual stations and in addition bilinear interpolation for multi-station data fills up spatial gaps also.The modified curve or gridded map is digitised with the values of VTEC including the interpolated values.The grid point hourly average VTEC values are given in tabulated form for real time applications.It is possible to observe graphical display of actual and model data for visualisation and decision making on the quality of VTEC model data at the control centre.Using elevations >50˚ and hourly time averages, over one station the satellites provide 3˚ × 3˚ latitude x longitude coverage of IPPs at a resolution of 1˚ × 1˚.The software has since been modified to provide VTEC model hourly average values also for two longitude zones covering the Indian region.The real-time model was submitted to GAGAN project and an evaluation was done by comparing the results with standard 5˚ × 5˚ grid point data received as part of satellite messages available with the project.The results of the author's real-time model for 13 August 2007 showed sharp gradients and pronounced EIA peak amplitudes compared to the smoothened GIVD values obtained through these messages.If the model values are also fitted with a smooth curve then the differences in GIVD values are reduced to ~25% -30%.However, further evaluation with more such comparisons would be needed.

Results
The GAGAN data for the period May, 2007-April, 2008 was provided in the format mentioned earlier in monthly CDs, each day of the month having one "MS Excel" file for each available station out of the 16 stations in GA-GAN network.

Single GPS Station Analysis
For one day's regional analysis the files could be selected based on number of stations required to be used for running the real-time TEC model.The resolution of the raw data provided is one minute and all satellite passes with elevation angles greater than 15˚ are used to derive the STEC and S 4 index values.For the purpose of getting better accuracies in the computation of VTEC, the data with satellite elevation angles greater than 50˚ are only used.This reduces the error due to the assumption in the mapping function with the whole ionosphere being condensed into a thin shell around 350 km altitude [22].The conversion from STEC to VTEC values is done only at sufficiently high elevation angles of greater than 50˚ to minimise the error due to this assumption of 350 km.Typical results of the computed VTECs on 15 May 2007 pertaining to Bangalore station are shown in Figure 3.
The panel 1) of the figure shows the latitude-longitude coverage of all satellite passes for the 24 hour period with selection of satellite elevation angles >50˚.Due to the application of this filter, the IPP's central latitudes and longitudes within the 3˚ × 3˚ grids are closer to the station coordinates as clear from this panel.The scintillation index S 4 is plotted at 1-min interval against IST (Indian Standard Time) starting at 5.30 hr.local time (0 hr. UT) as shown in panel 2).Any abnormal variations of the scintillation index serves as an alert to check for any impacts on the TEC values.This is a good cross-check to evaluate unusual changes in VTEC values.Panel 3) plots the 1-min VTEC values from different satellite passes and shows the cubic spline fitted mean curve using all values between 12 -14˚N latitude and 77 -79˚E longitude.Based on the spline fitting the hourly average values (like average values between 5.30 hr -6.30 hr. is assigned at the 6.00 hr.value) of VTEC are plotted in panel 4) corresponding to the 3˚ × 3˚ grid average values over Bangalore.The same values of VTEC are plotted as pixels with latitudinal variation (for all longitude values between 77 -79˚E) without interpolation in panel 5) and with cubic spline interpolation (for time gaps) shown in panel 6).

Multi-Station Analysis and Latitudinal Variation of VTEC
For running the multi-station real-time model the observations (individual data files of available network stations) are run one station after another so that the VTECs can be estimated for the previous 24 hours.Each stations geographic coordinates are given as input values along with the selection of elevation angle which is recommended to be kept at >50˚.Time averaging is done on hourly basis as a default setting.Each stations output is provided with tabulated hourly average data and graphical presentation of results as shown already in Figure 3.After completing all the individual stations data analysis the computer programme continues to carry out integrated analysis by taking into account the 3˚ × 3˚ grid values at 1˚ × 1˚ resolution (of each station) and time versus latitude pixels at 1˚ latitudinal resolution are generated as output tabulated data along with graphical presentation of pixelated VTEC values for instant monitoring separately for two longitude zones, i.e., 73 -83˚E and 83 -93˚E as well as for the entire longitude range of 73 -93˚E.
Figure 4 shows a sample of such a result from a multistation data analysis on 13 August 2007 using only the best data available for 9 stations fairly well distributed in latitude and longitude over the Indian region.If the data of Trivandrum (8.5˚N, 76.9˚E) station is also available then the model output values would cover up to 6˚N lowest latitude instead of 11˚N.The colour code has a resolution of ~5 TEC units.The tabulated data gives the actual VTEC values with 1 TEC unit resolution.The white pixels indicate data gaps not possible to be interpolated by the present algorithm.Even taking the 1˚ × 1˚ resolution in the 3˚ × 3˚ coverage from individual stations some latitude-longitude pixels may still not be filled and this situation would improve with a denser and more uniform distribution of GPS receiver stations.Due to the peculiarities of the landmass distribution of the region and inadequate uniform distribution of stations particularly with longitude, there is a good option to consider the region as one longitude zone and concentrate only on latitudinal variation at the cost of some loss of accuracy due to possible longitudinal variations in VTEC discussed later.

Day-to-Day Variability of VTEC
The real-time model is further evaluated in two steps.In the first step the VTEC contours for two consecutive days of 13 and 14 August 2007 are plotted in Figure 5.
These contours are derived from the pixelated data of the type shown in Figure 4 with additional bilinear interpolation to fill up some spatial gaps.The values of the VTEC for the two days are similar but for a small component of day to day variability even when the F10.7 solar flux values remained the same at 65 units for both days.This demonstrates that the model is sensitive to bring out non solar activity related small day to day variations.So the model can be easily adopted to study the day to day variability of VTEC over Indian region using such observations during different seasons and years to separate solar and non-solar activity related causes of such variations.Figure 5 also shows that while the low latitude ionospheric anomaly peaks clearly appear on both days there are differences in details on peak latitude(s) and their time(s) of occurrence.

Longitudinal Variation of VTEC
In the second step consideration is given to distribution of the GAGAN TEC stations over a longitude span of 20˚.This requires that the model results are checked with respect to any longitudinal variations.Since longitudinal variations are normally expected to be small, the test is carried out by dividing the longitudes into two zones, one between 73 -83˚E and the other between 83 -93˚E.The 0530 hr.IST (equivalent to 0 hr.UT) is based on the 82.5˚E longitude time which is at the central point of the above longitude span.The real-time model automatically provides for VTEC maps of the same day for these two longitude zones for comparison.Such time-latitude contour plots for two longitude zones are shown for 13 and 14 August 2007 in Figure 6.The model automatically selects the stations for the two longitude groups out of the available inputs.In the example for comparison of the two contours the latitude coverage is restricted between 17 -28˚N.While on 13 th August there is some indication of slight time shift of early occurrence in the longitude zone 83 -93˚E, the changes are only marginal and hence using the VTEC values irrespective of longitude differences within 73 -93˚E may serve the purpose of a model of this region.However, for GIDV values a suitable decision can be taken at the control/monitor whether to split the VTEC values into 2 longitude zones or not based on real time situation and results of previous model runs every half an hour.The present model provides for selection of this option as the contours for the total longitude range are also presented along with the other two plots.The longitudinal anomalies may be interesting to study during solar energetic events or space weather related changes which propagate with earth's rotation affecting different longitude zones at different times.

Sub-Latitudinal Variation of TEC
Considering that longitudinal variations of VTEC within the GPS station network are not very pronounced, the  of VTEC during August are small occurring much after the sunset, these are very pronounced in June occurring near the sunset.To further quantify the relative amplitudes of the VTEC values the results of the same pairs of days are plotted individual day wise and shown in Figure 8.It can be noted that in all the four days the amplitudes of VTEC peaks are strongest for the anomaly latitude group followed by the equatorial and low latitude groups.The overall values of VTEC are marginally higher for June as compared to August.Hence it is possible to work out a monthly average pattern or shape of the diurnal curve of VTEC which could serve as very useful reference for the daily predictable model outputs.

Conclusions
1) Ionospheric research has been vigorously pursued over the Indian region in well-coordinated manner since the IGY period.Apart from investigating various scientific problems peculiar to the low latitude related variations covering the equatorial electrojet and EIA regions, the studies helped in direct applications of the results to HF radio communication in earlier days.Through a fast pace of development in the area of upper atmosphere and space science and technological growth of microwave and satellite communications, the emphasis has now shifted to develop reliable ionospheric models with wide ranging applications in satellite communication and navigation systems.
2) The requirement of global and regional ionospheric models led to the new observational concept of networking of ground stations for near simultaneous observations needed for high accuracy in e.g.global positioning and radar tracking applications.Luckily, the opportunity for such multi station observations is provided by the high technology system itself like in the case of measuring TEC by dual frequency GPS receivers.
3) With the advent of GNSS system and its applications to precise location determination, e.g. for aircraft landing etc., many countries initiated setting up ground based receiving stations at regional scales with a view to check available ionospheric TEC models as well as improve the accuracies based on new observations.India initiated the GAGAN project for the purpose with one of the objectives to develop a suitable real-time regional model to estimate range errors over the Indian low latitude region.The need for such a project was highlighted in view of the fact that available global models are not adequate to cater to the highly variable TEC values with sharp gradients peculiar to this region.
4) The author has developed a real-time VTEC model based on the data provided from a network of ground stations by the GAGAN project.For each run of the programme written in MATLAB, the algorithm uses the 24 hr.data with one minute resolution from each of the operational TEC stations spread over in latitudes (8.5 -31˚N) and longitudes (73 -93˚E), selects only those GPS satellite passes which have elevation angles >50˚, uses a 6) While the real-time model has gone through many modifications and revisions, the variants of the basic core programme provide avenues for very detailed scientific investigations with the new sets of homogenised data collected under the GAGAN project.As an illustration, the results are presented by dividing the latitude range into 3 sub-groups i.e., equatorial (10 -17˚N), anomaly (17 -25˚N) and low latitude (25 -32˚N).Apart from the expected peaks of mean VTEC values for the anomaly latitude group compared to the other two, there is a marked change in the overall diurnal pattern of VTEC variation or the shape of the curve between August and June months.Such analysis for all the months would provide the trends in VTEC variations which would be useful for short term predictions.In addition, the impact of variations due to different solar activity periods may be studied using the present model with additional data.Such scientific results could provide deep insights to the variability of TEC over the Indian region and the information can be fed back for improving the real-time model with new elements of short term prediction features to reduce the dependence on large number of ground observation stations.

Figure 2 .
Figure 2. Distribution of dual frequency GPS ground receiver/TEC stations under GAGAN project.

Figure 3 .
Figure 3. Results of observed and interpolated VTEC values with related GPS satellite parameters for one GAGAN network station (Bangalore) on 15 May 2007.The elevation angles of GPS satellites chosen is >50˚.The panels (a) to (c) are based on 1-min raw data and (d) to (f) are analysed results after subjecting to cubic spline interpolation/fitting and hourly averages.

Figure 4 .
Figure 4. Multi station analysis using the real-time TEC model using only 9 stations data to generate VTEC map at 1˚ latitude resolution for any longitude between 73 -93˚E covering the Indian region.The colour resolution is for ~5 TEC units and the white pixels denote no data.Only the cubic spline interpolation is used to cover time gaps and spatial interpolation is restricted only at individual station pixel level to achieve 1˚ resolution within 3˚ latitude values.

Figure 5 .
Figure 5. Time-Latitude contour plots of VTEC applicable to the whole longitude range (73 -93˚E) of the GAGAN network stations using the real-time model for 13 th and 14 th August 2007.The white areas represent no results.

Figure 7 .
Figure 7.Diurnal variation of VTEC on few selected days divided into 3 sub-latitude regions of equatorial (10 -17˚N), anomaly (17 -25˚N) and low latitude (25 -32˚N) all for the whole longitude range (73 -93˚E).The days are taken in pairs to show the effect of day to day variability and changes in the shapes of the curves.

Figure 8 .
Figure 8.Diurnal variation of VTEC for the 3 latitude groups as in Figure 7 for individual days in summer months.mapping function to convert the hourly mean slant STEC to VTEC values provided at hourly intervals and for three longitude blocks i.e., for the entire longitude range 73 -93˚E and for 73 -83˚E and 83 -93˚E.The results are given in output tables at 1˚ latitudinal resolution as well as presented graphically as time-latitude maps and contours.Cubic spline interpolation technique is mainly used for filling up results in time gaps and bilinear interpolation for spatial gaps.5) Sample results of few days in May, June and August of 2007 are presented which could be available to the central monitoring station in real time.There is good agreement in VTEC values for two sample pairs of closeby days.The VTEC variations with the change in longitude blocks are very small at least for these selected days of very low sunspot activity.6)While the real-time model has gone through many modifications and revisions, the variants of the basic core programme provide avenues for very detailed scientific investigations with the new sets of homogenised data collected under the GAGAN project.As an illustration, the results are presented by dividing the latitude range into 3 sub-groups i.e., equatorial (10 -17˚N), anomaly (17 -25˚N) and low latitude (25 -32˚N).Apart from the expected peaks of mean VTEC values for the anomaly latitude group compared to the other two, there is a marked change in the overall diurnal pattern of VTEC variation or the shape of the curve between August and June months.Such analysis for all the months would provide the trends in VTEC variations which would be useful for