Evaluation of a global model of ionospheric slab thickness for foF2 estimation during geomagnetic storm

This study aims to evaluate the performance of the global model of ionospheric slab thickness (GMIST) in terms of F2 layer critical frequency (foF2) estimation during geomagnetic disturbed conditions. Hourly values of foF2 as obtained from ionosonde stations located at equatorial, low- and mid-latitude regions are compared with the corresponding GMIST and IRI-STORM modeled values. For this purpose, the correlation coefficient, daily mean, root mean square error and improvement percentage are calculated at different regions and geomagnetic disturbance levels. The results show that GMIST is more accurate than IRI-STORM model in terms of foF2 estimation at low- and mid-latitude regions, while at equatorial areas GMIST is less accurate during geomagnetically disturbed and quiet conditions.


Introduction
The total electron content (TEC), measured by means of navigation satellite systems such as GPS, GLONASS, GALILEO, has turned into a key parameter in the characterization and monitoring of ionosphere on a global scale.However from the point of view of ionospheric radio propagation the most important characteristic is the critical frequency foF2, whose value directly defines an optimum working frequency for near veritcal incidence skywave (NVIS) propagation.The ionospheric density profile and its maximum electron density (NmF2) is directly related to foF2 which is be regularly measured by ionosondes.The spatial coverage of foF2 ionosonde observations is not globally uniform due to the high cost of acquiring and maintaining ionosondes.This limitation can be surpassed considering the abundance of GPS receiver networks that provide TEC observations which can be used jointly with ionospheric slab thickness models to improve the foF2 predictability and therefore to provide support for various HF communication systems, radio amateurs and broadcaster operators.Houminer [1997] has reviewed the use of GPS data in relation to short-term foF2 prediction enhancement using one ionosonde station located in Cyprus and concluded that the high correlation between TEC GPS with foF2 (varying from 0.53 to 0.78) allows the use of GPS data for near real time foF2 map updating, both in temporal (short-term, medium and long-term) and in spatial (local, regional and global) terms.Kouris [2004] concluded that there is a high correlation between the TEC and the square of foF2 (around 0.8), although ionospheric hysteresis is detected before and after noon time.Ma and Maruyama [2002], based on GPS TEC data and Kokubunji ionosonde foF2 data, reported that seasonal variations of ionospheric slab thickness during a daytime exhibits opposite phase with respect to the nighttime variation.Ionospheric slab thickness depends on three parameters; latitude, season, and solar activity.In the summer, ionospheric slab thickness is greater during daytime than in the night.The ionospheric slab thickness depends linearly on the solar activity.Jayacandran et al. [2004] concluded that the diurnal variation of the ionospheric slab thickness during low solar activity period is generally characterized by higher value in the night than during the day, while the situation is reversed during high solar activity at the non-auroral region.In auroral areas, the nighttime slab thickness is higher than the daytime one during both low and high solar activity periods.Variability of ionospheric slab thickness is higher in the night than during the day in the two phases of solar activity for mid-, low-and high-latitude ionosphere.Suvannasang et al. [2008] con-cluded that the diurnal variation of the ionospheric slab thickness reaches its maximum value at 05:00 LT.The average daytime slab thickness is minimized at around equinox, whereas the average nighttime slab thickness is maximized during summer.Gerzen et al. [2013] have used Neustrelitz global TEC model (NTCM-GL) [Jakowski et al. 2011] to simulate TEC and the Neustrelitz Peak Density Model (NPDM) [Hoque and Jakowski 2011] to generate NmF2.Based on these two ionospheric parameters, they model ionospheric slab thickness.Then, ionospheric slab thickness model was employed for NmF2 estimation from operational GNSS global TEC model.Validation of the NmF2 reconstruction from GNSS TEC model shows that both the model and the reconstruction approach produce improved results, however, the reconstruction has a factor of 1.75 better mean value during low to middle solar activity period.Still, NmF2 reconstruction has not been validated during an ionospheric storm.Muslim [2010] has developed a global model of ionospheric slab thickness (GMIST) by using global ionospheric map (GIM) model and measured foF2 ionosonde values.The aim of this study is to evaluate the GMIST model performance in terms of foF2 estimation during disturbed geomagnetic conditions.In particular, the capacity of the model to estimate foF2 during six geomagnetic storm events (moderate, strong, severe, and great) during the period 2000-2014 was assessed.

Ionospheric slab thickness
TEC can be expressed as an integral of electron density ranging from a lower limit to an upper limit of the ionospheric layer.Considering electrons from 60 km to 20,000 km TEC can be written as If the single layer model of the ionosphere is assumed at an altitude of 350 km and NmF2 is the maximum electron density of the ionosphere, Equation (1) can be expressed as where x is equivalent to the ionospheric slab thickness.
Ionospheric characteristics are extracted from ionograms that are inverted to produce altitude profiles of electron density.The critical frequency of the F2 layer (foF2) is related to the maximum electron density in the F2 layer NmF2 as follows: So Equation (2) becomes and ionospheric slab thickness can be expressed in term of foF2 by where TEC is measured in TECU (express this in electrons per cubic meter), foF2 in MHz, and x in km.

Global model of ionospheric slab thickness
GMIST formulations are based on the model of simplified low-latitude region of ionosphere MSILRI [Muslim et al. 2007, Muslim 2010], which are harmonic functions that describe ionospheric diurnal variation, and use polynomial function for representing the latitudinal variation of the ionosphere, and a linear function to capture the ionospheric slab thickness response to the solar activity.Since the thickness model is a global model that includes high, mid, low and equatorial latitudes, the geographical latitude dependence of MSILRI is replaced by geomagnetic latitude dependence.
Considering data availability for 2 days ahead prediction of foF2, Muslim [2010] has used TEC data derived from global ionospheric maps that are produced every day at the Center for Orbit Determination in Europe (CODE GIM).CODE GIM based on GPS data as obtained from about of 150 IGS receivers and other institutions.TEC is modeled in a solar geomagnetic reference frame using a spherical harmonic expansion up to degree and order of 15.GIM time resolution is 2 hours and spatial resolution is 2.5° and 5° of latitude and longitude, respectively.GIM models can be obtained from ftp://.ftp.unibe.ch/aiub/CODE.In order to develop the ionospheric slab thickness model, GIM data during the period 2000-2009 were employed.The ionospheric foF2 data were obtained from ionosonde observations located in Sumedang, Darwin, Vanimo, Singapore, Manila and other locations as listed in Table 1 and shown in Figure 2.
The GMIST model is developed considering the following four assumptions: (a) linear response of foF2, TEC and slab thickness to solar activity for both descending and ascending solar cycle, (b) no longitudinal dependence of the three ionospheric parameter at same local time, (c) conversion of geographical to geomagnetic coordinates does not affect the local time variation of ionosphere and (d) slab thickness is constant during variation of foF2 and TEC for both quiet and disturbed geomagnetic conditions.The last assumption allows for slab thickness model application for short term foF2 prediction using real-time TEC data.
GMIST employs the MSILRI formulation so as to .
. .describe the ionospheric variation with latitude using geomagnetic latitude instead of geographical latitudinal representation.At a certain ionosonde station, TEC and foF2 response to solar activity index R12 are expressed as R12 is used in GMIST model due to the predictibility and availability of long historical data that make possible the recontruction of past time foF2 data from TEC and vice-versa.Using Equations ( 6) and ( 7) at a month for all local time, the ionospheric slab thickness at station x of Equation ( 5) is written as where x refers the station and t refers to local time.
Eleven order harmonic analysis of the estimated slab thickness from t =0 until 23, results to 23 coefficiens of harmonic function Harmonic analysis of slab thickness for all stations results to X sets of harmonic coefficients Cn and Dn which can be represented as a function of geomagnetic latitude using From Equations ( 10) and ( 11) we have two sets of polynomial coefficients of c n for R12 = 0 and 100 as well as polynomial coefficients of s n .Using the two sets of the polynomial coefficients, the linear regression of polynomial coefficients to solar activity is expressed by equations GMIST procedures adopt the MSILRI procedures.MSILRI procedures adopt the simplified ionospheric re- ( (13) gional model (SIRM) procedures [Zolesi et al. 1996].
Figure 1 shows the GMIST procedures as follows: (a) From a set of monthly mean TEC at local time t, and monthly smoothed R12, the linear regression of TEC against R12 results to coefficients A0 and A1 using Equation ( 6).Linear regression of monthly median foF2 versus R12 results to coefficients B0 and B1 with Equation ( 7).
(b) Using linear regression coefficients, TEC-mod100 and foF2mod100 arecalculated for R12 = 100, that represent TEC and foF2 model for high solar activity as well as TECmod0 and foF2mod0 for low solar activity R12 = 0.
(c) Then the ionospheric slab thickness at both solar activities level, x0 and x100 is obtained using Equation ( 9).
(d) Repeating procedures (a) until (c) for all local time from 0 until 23 LT, the 24 values of ionospheric slab thickness at R12=0 and R12=100 are produced.Using harmonic analysis of the 11th order, 23 coefficients of the harmonic function are obtained using Equation ( 9).
(e) Repeating procedures (d) for all stations, each harmonic coefficients can be represented by using the polynomial function with the polynomial regression for R12=0 and 100.
(f ) Polynomial coefficients of harmonics coefficients of Equations ( 10) and ( 11) can be represented as R12 function using the linear response assumption of the ionosphere to solar activity.The last procedure results to coefficients a and b.

Data and methodology
In order to evaluate GMIST capability to produce foF2 values during disturbed geomagnetic conditions, we have used foF2 data as obtained from ionosonde sta-  Table 2. List of geomagnetic storm events used in GMIST evaluation.
tions located at mid, low, high latitude and equatorial areas (see Table 3 and Figure 2) during six geomagnetic storm events that are shown in Table 2.These storm events were classified as moderate, strong, severe, and great with the aid of Dst geomagnetic index [Loewe and Prölss 1997].As reference days characterized by quiet ionospheric conditions, we considered one or two days before the geomagnetic storm event.The critical frequency of F2 layer was calculated using GMIST model as derived from CODE GIM.The accuracy of the GMIST foF2 has been quantified by calculating the daily RMSE between the model and observed foF2 values and by comparing it with the foF2 RMSE value calculated by the IRI-STORM model The equations of RMSE and percentage improvement that were used are:  ( where m i =the model output (GMIST or IRI), o i =the observed values, I=the percentage improvement The RMSE and percentage improvement, I, are calculated for different geomagnetic latitude and geomagnetic disturbance levels.

Results
Figure 3 shows the diurnal variability of Dst index and estimated foF2 by GMIST, IRI and observations data for the quiet day (Dst > −30 nT) of July 14, 2000, and the severe geomagnetic storm (Dst < −300 nT) of July 15-16, 2000.The vertical blue line in Figure 3b represents the time of geomagnetic storm peak.During geomagnetic storm, the GMIST model underestimates the foF2 observations and the IRI foF2 modeled values especially the foF2 minimum occurring before sunrise at the lowlatitude ionosphere (Figure 3b, right-lower panel) and mid-latitude Southern Hemisphere (Figure 3b, left-lower panel).The largest negative ionospheric storm occurred on July 15, 2000, following the peak of geomagnetic storm time.GMIST modeled foF2 values shows a good agreement with observational foF2 values during quiet and geomagnetic storm in the Northern Hemisphere mid-latitude region (Figure 3b, upper panel).
In Figure 4, the scatter diagrams between modeled (GMIST and IRI) and observed ionosonde foF2 values are presented.In addition, the correlation coefficient (R), the foF2 mean (M) and the RMSE values for a quiet day ( July 14, 2000) and during the geomagnetic storm of July 15-16, 2000, are shown in Figures 3a and 3b respectively.As it can be observed, the correlation coefficient between GMIST and foF2 observations are high (R > 0.9) both during quiet and disturbed conditions, in all areas except for the low-latitude areas (Figure 4b, right-lower panel).However, during geomagnetic storm, the mean values of GMIST and IRI foF2 are lower at low latitude and higher at mid latitude than observations (Figure 4b, right-upper panel).As it can be seen from Figure 4 (right-lower panel), the mean foF2 modeled value at low latitude is around 9.59 MHz 4, or 0.9 MHz lower than the foF2 observed value.Contrary to low-latitude region, mean foF2 modeled values at midlatitude are higher than observed foF2 values, with modeled foF2 being equal to 6.18 MHz whereas mean observed foF2 values being around 5.7 MHz (Figure 4b, right-upper panel).Compared to the mean values of foF2 during quiet day (Figure 4a), the GMIST foF2 mean values during geomagnetic storm (Figure 4b) are about 0.4 MHz higher at mid-latitude areas and 0.2 MHz lower at low-latitude areas.
The Dst index variation, and foF2 values as generated from GMIST and IRI models and observed from ionosondes during the great geomagnetic storm of October 29-31, 2003, are presented in Figure 5.As it can be seen, the GMIST foF2 is in good agreement with observations for the low-latitude region (Figure 5b, rightlower panel) especially during positive ionospheric storm, while after the peak of the storm, the GMSIT model overestimates foF2 observed values.In the Southern Hemisphere, negative ionospheric storms occurred at low-and mid-latitude regions (Figure 5b  GMIST corresponds to the equatorial region both for quiet days and geomagnetic storm (Figures 6a and 6b, right-upper panels).High overestimations of GMIST foF2 values were detected in all stations during this severe geomagnetic storm.
In Figure 7, the modeled (using GMIST and IRI-STORM) as well as the observed foF2 variations during the moderate geomagnetic storm of April 4-6, 2010, are demonstrated.Qualitatively, GMIST performance is superior to IRI STORM model at low-and mid-latitude regions, while it is less accurate than IRI STORM model at the equatorial region.The largest RMSE is also detected in the equatorial region (Figure 7c, rightlower panel).
During the geomagnetic storm of July 15-16, 2000, overall, the GMIST and IRI RMSE are increased.The RMSE between modeled (GMIST and IRI) observations for selected geomagnetic storm events and selected ionosonde stations are showed in Table 4.The improvements of GMIST foF2 over IRI STORM foF2 values were also calculated.In the same table, the RMSE and percentage improvement were grouped into 3 ionospheric regions; mid-latitude, low-latitude and equatorial regions and separated between quiet and storm time.As it can be seen, in general GMIST foF2 improvement is positive (efficient) for mid-and low-latitude areas but it has negative value (or inefficient) at the equatorial region during both quiet and disturbed geomagnetic conditions.Significant improvement is noted at low-latitude regions where the GMIST improvement accuracy over IRI STORM foF2 is by 29.73% during geomagnetic storm and by 44.44% for quiet days.At the mid-latitude region, the GMIST improvement efficiency is significant only during geomagnetic storm where percentage improvement is approximately 14.57% while during quiet day the percentage improvement is minor (5.27%).
More analytically, GMIST model demonstrates higher improvement at the mid-latitude regions, in all stations apart from Learmonth during the storm event of July 15-16, 2000, with the negative percentage improvement value −51.52% approximately.The GMIST estimate foF2 over Learmonth is lower than observational data.This was observed during the severe negative ionospheric storm event of July 15-16, 2000 (max Dst −301 nT).
Concerning Nicosia station, during the recent moderate geomagnetic storm of April 12, 2014, the GMIST model RMSE is around 0.1 MHz whereas the IRI STORM model RMSE is about 0.09 MHz, and correlation coefficient between modeled and observations val-

| -a
| -b ues is higher than 0.95 both for GMIST and IRI-STORM model.The corresponding percentage improvement has a negative value, indicating that IRI-STORM model can approximate well the observations.On the contrary, during the geomagnetic storm that occurred on March 7-9, 2012, GMIST model shows a positive foF2 improvement around 13%. Overall, at Nicosia station the GMIST model performance in terms of foF2 calculation is improved by 3%.On the opposite, the GMIST for Chilton during the October 2003 great geomagnetic storm was not accurate as the respective percentage improvement was negative and around −13.63%.However, during the moderate geomagnetic storm GMIST model was capable of reproducing foF2 more accurately than IRI STORM model at about 33.38%.
As regards the foF2 estimation at low latitudes a significant improvement of GMIST model as compared with the IRI STORM model was concluded.Apart from the geomagnetic storm of July 2000 where a negative improvement (−6.95 %) was detected, GMIST model improvement during geomagnetic storms (Table 2) is high (>23%).In general, the GMIST model can perform better than the IRI-STORM model at the low-latitude regions during quiet geomagnetic conditions with the percentage improvement reaching 44.44%.Opposite to low-latitude regions, at equatorial regions, GMIST performance in terms of foF2 estimation is worse compared to the IRI STORM model performance during both storm and quiet day as seen in Table 4.This ineffectiveness is caused by the overestimation of GMIST over IRI STORM and observations as shown in right-upper panel of Figure 4b.During the negative In Table 4, three to fourth negative percentage improvement at mid and low latitude are occurred during the severe geomagnetic storm.Figure 8a shows that the GMIST RMSE from observational data is increase with increasing geomagnetic disturbance level at mid, low latitude, and the equatorial region.The bottom panel of Figure 7b shows that the GMIST RMSE is lowest at mid-latitude at about 42.12°ML.Largest RMSE is at low latitude for both models, followed by RMSE at the high/mid-latitude region of about 53.81°ML.

Discussions
We found that the GMIST foF2 estimation improvement over IRI STORM during the geomagnetic storm is significant for the low-latitude region, moderate effective for middle latitude and ineffective for equatorial region.This is caused by latitudinal variability of ionospheric response to the geomagnetic storm.High/mid-latitude, low-latitude, and equatorial ionosphere are affected significantly by the geomagnetic storm that caused the slab thickness changing more than at mid latitude.Constant slab thickness assumption of GMIST caused foF2 estimation is also a  contributing factor for errors.
The GMIST foF2 improvement is highest during geomagnetic storm due to the use of low-latitude Sumedang ionosonde data for GMIST modeling.However, GMIST RMSE over Tucuman low-latitude region that is not included for GMIST modeling is highest compared to mid-latitude and equatorial ionosphere.
Ineffectiveness of GMIST foF2 for equatorial regions can be caused by data used for GMIST modeling that only include the Asian sector for equatorial region while validation uses equatorial data at the American sector.This is caused by large longitudinal variability between Asian and American sector.GMIST also use assumptions of no longitudinal dependence on ionospheric parameters in a local time reference frame.Ignor-ing the longitudinal variability of GMIST contributes significantly to the ineffectiveness of GMIST foF2 not only for an equatorial region but also for mid-and low-latitude region.Ignoring longitudinal variability also caused overestimation of GMIST during quiet days.
Compared to geomagnetic storm condition, GMIST improvement over IRI STORM for quiet days is more effective for low-latitude region.Although less effective for mid-latitude region the GMIST improvement is 5.27% over IRI STORM.
Based on model performance, the GMIST can be used as an alternative model for foF2 prediction during quiet and storm geomagnetic condition at mid-and low-latitude region.However, significant revision is required to improve GMIST accuracy.Revisions of GMIST will include longitudinal variability and geomagnetic disturbance index as an input parameter.

Conclusions
Evaluation of GMIST foF2 estimation was carried out using ionosonde data for 5 geomagnetic storm events ranging from moderate to great.The results show that for the mid-latitude region the GMIST foF2 improves a storm time foF2 estimation over IRI STORM by 14.57%, and for low latitude by 29.73%.For the equatorial region, GMIST foF2 improvement over IRI STORM is insignificant due to GMIST estimating foF2 higher than observations even during quiet days.
To improve the accuracy of the GMIST foF2 estimation for equatorial regions, more foF2 data over the equatorare needed, and longitudinal variability should also be examined.To improve GMIST accuracy during severe and great geomagnetic storms, the modeling procedures should consider geomagnetic disturbances index as a paremeter input to GMIST.

Figure 2 .
Figure 2. Map of ionospheric stations used for GMIST development and evaluation.

Figure 3 .
Figure 3.Diurnal variation of (a) Dst index and (b) foF2 estimated by GMIST, IRI and ionosonde observations on July 15-16, 2000, at ionospheric stations of different regions of the globe.The geomagnetic latitude and geographic longitude of each station is given at the low-right edge of each plot.
, lower panel).Close to equator, GMIST model underestimates foF2 observations around the minimum of diurnal variation on October 30, 2003, while it overestimates foF2 values near the diurnal foF2 maximum (Figure 5b, right-upper panel).In general, GMIST overestimates foF2 at high/mid-latitude areas as seen in the left-upper panel of Figure 5b.The scatter diagram between modeled (GMIST and IRI) and observed foF2 values are shown in Figure 6.The highest RMSE for IONOSPHERIC SLAB THICKNESS FOR foF2

Figure 4 .
Figure 4. Scatter diagram between modeled (GMIST and IRI STORM) and observed foF2 values during (a) quiet day and (b) geomagnetic storm on July 15-16, 2000, at ionospheric stations of different regions of the globe.The geomagnetic latitude and geographic longitude of each station is given at the upper-right edge of each plot.

Figure 5 .
Figure 5. (a) Dst index during period of great geomagnetic storm and (b) estimated foF2 from GMIST, IRI, and ionosonde observations at ionospheric stations of different regions of the globe.The geomagnetic latitude and geographic longitude of each station is given at the lowright edge of each plot.

Figure 6 .
Figure 6.Scatter diagram of estimated foF2 by GMIST and IRI STORM from observations at ionospheric stations of different regions of the globe.The geomagnetic latitude and geographic longitude of each station is given at the low-right edge of each plot.

Figure 7 .
Figure 7. (a) Dst index, (b) comparison of GMIST foF2, IRI model and observations during April 4-6, 2010, and (c) scatter diagram GMIST foF2 and IRI vs foF2 data at ionospheric stations of different regions of the globe.The geomagnetic latitude and geographic longitude of each station is given at the low-right edge of each plot.

Figure 8 .
Figure 8.Effect of the maximum deviation of Dst index on RMSE for mid-, low-latitude and equatorial region (upper panel) and geomagnetic latitude variation of RMSE (bottom panel).

Table 1 .
List of ionosonde stations used for ionospheric slab thickness model development.

Table 3 .
List of ionosonde stations used for the evaluation of GMIST.

Table 4 .
RMSE and Improvement GMIST and IRI during storm and quiet days.