A new global model for the ionospheric F 2 peak height for radio wave propagation

The F2-layer peak density heightmF2 is one of the most important ionospheric parameters characterizing HF propagation conditions. Therefore, the ability to model and predict the spatial and temporal variations of the peak electron density height is of great use for both ionospheric research and radio frequency planning and operation. For global hmF2 modelling we present a nonlinear model approach with 13 model coefficients and a few empirically fixed parameters. The model approach describes the temporal and spatial dependencies of hmF2 on global scale. For determining the 13 model coefficients, we apply this model approach to a large quantity of global hmF2 observational data obtained from GNSS radio occultation measurements onboard CHAMP, GRACE and COSMIC satellites and data from 69 worldwide ionosonde stations. We have found that the model fits to these input data with the same root mean squared (RMS) and standard deviations of 10 %. In comparison with the electron density NeQuick model, the proposed Neustrelitz globalhmF2 model (Neustrelitz Peak Height Model – NPHM) shows percentage RMS deviations of about 13 % and 12 % from the observational data during high and low solar activity conditions, respectively, whereas the corresponding deviations for the NeQuick model are found 18 % and 16 %, respectively.


Introduction
The ionospheric F2-layer is primarily responsible for the reflection of high frequency (HF) radio waves in the ionosphere.Thus, the F2-layer peak density height hmF2 is one of the most important parameters needed for radio frequency planning and spectrum management.The regular variation of the solar radiation with the solar zenith angle causes temporal and spatial variations of hmF2.Depending on the solar activity, daytime and season, the peak density height may range from 350 to 500 km at equatorial latitudes and from 250 to 350 km at mid-latitudes.Additionally, there is a strong dependency of hmF2 on dynamic forces such as electric fields and neutral winds.Due to regular and irregular variations of the bottomside plasma density closely related to the NmF2 and hmF2 variations, the terrestrial signal transmission may be interrupted or even lost; moreover, the transmission coverage may be affected due to up or down lifting of the ionospheric plasma, changing the hmF2 height.
Furthermore, the Earth-space transionospheric communication can also benefit from the knowledge of the hmF2 and NmF2.As an example, GNSS (global navigation satellite system) positioning can be improved by mitigating higher order ionospheric propagation effects such as ray path bending errors using NmF2 and hmF2 information (see Hoque andJakowski, 2008, 2011a).Again, since F2-layer peak is a key anchor point to construct ionospheric electron density profiles, NmF2 and corresponding hmF2 are the most important parameters in empirical ionospheric modelling.The accuracy of the peak height is crucial in some other applications too, such as inferring the neutral wind (e.g.Zhang et al., 2003).
To develop hmF2 prediction model, early work was done by Shimazaki (1955), Bradley and Dudeney (1973), Bilitza et al. (1979), and Dudeney (1983) utilizing the propagation factor M(3000)F2.The M(3000)F2 can be deduced from vertical-incidence ionograms using standard methods (Piggott andRawer, 1972, 1978).The propagation factor M(3000)F2 is related to the maximum usable frequency MUF(3000) by M(3000)F2 = MUF(3000)/foF2, where MUF(3000) is defined as the highest frequency at which a radio wave can be received over a distance of Published by Copernicus Publications on behalf of the European Geosciences Union.3000 km after reflection in the ionosphere (Bradley and Dudeney, 1973).Shimazaki (1955) found that hmF2 is inversely related to M(3000)F2.Bradley and Dudeney (1973), Bilitza et al. (1979), and Dudeney (1983) obtained better results after considering dependency of hmF2 on the ratio of critical frequencies in the F2 and E layer foF2/foE, sunspot number and geomagnetic latitude, in addition to its dependency on the M(3000)F2.In another approach, McNamara et al. (1987), Kishcha and Kochenova (1996) computed hmF2 directly from ionosonde measurements using true height analysis.
Similar work was done by Jones and Obitts (1970), Rush et al. (1983Rush et al. ( , 1984)), Fox and McNamara (1988), Bilitza et al. (1990), Bilitza (2001), and Bilitza and Reinisch (2008) in developing global ionospheric parameter models.Considering the importance of ionospheric characteristics in radio frequency planning, the International Telecommunications Union (ITU) issued a standard set of ionospheric parameter models (CCIR, 1967;ITU-R, 1997) on advice from its International Radio Consultative Committee (CCIR), presently from its Radio-Communication Sector (ITU-R).The CCIR models and related software are available via ITU.The CCIR model consists of 24 maps, each one containing 441 coefficients for one month of the year and one of the two levels of solar activity, R12 = 10 and 100, where R12 is the 12-month running mean of the monthly sunspot number (1764 coefficients in all) (Jones andGallet, 1962, 1965).Many empirical models such as the International Reference Ionosphere (IRI) model (Bilitza and Reinisch, 2008), the NeQuick model (Radicella and Leitinger, 2001;Coisson et al., 2006;Nava et al., 2008) and even some theoretical models use the CCIR maps for foF2 and hmF2 estimations.The IRI model uses modified Bilitza et al. (1979) equations for hmF2 estimation in any place and time using CCIR maps.
Another approach of ionospheric parameter modelling based on empirical orthogonal function (EOF) analysis of the observational data set is studied by Liu et al. (2008), and Zhang et al. (2009). Furthermore, Gulyaeva et al. (2008) recently developed a numerical model of hmF2 using about 90 000 electron density profiles derived from observations taken by topside sounder satellites ISIS1, ISIS2, IK19 and Cosmos-1809during 1969-1987.In this paper, we derive an empirical hmF2 model based on non-linear least squares technique.For this, we considered a set of nonlinear equations with 13 polynomial coefficients describing the regular variations of the peak density height.The polynomial coefficients are derived from a nonlinear fit with hmF2 measurements in least squares sense.For this, we used two types of measurement data, namely space-based GNSS ionospheric radio occultation (IRO) measurements and ground-based ionosonde observations.CHAMP (CHAllenging Minisatellite Payload) IRO data are used for inputs under both solar maximum and minimum conditions (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008), whereas GRACE (Gravity Recovery And Climate Experiment) and COSMIC (Constellation Observing System for Meteorology, Ionosphere and Climate, also known as FORMOSAT-3) IRO data are used for the low solar activity period 2006-2010.Additionally, we used a large database of propagation factors M(3000)F2 and hmF2, collected through a worldwide network of 69 ground ionosondes over the last 60 years.

Data sources
The most powerful source of hmF2 data used in the present study is the IRO observation.We used hmF2 estimates reconstructed from IRO measurements onboard CHAMP, GRACE and COSMIC satellite missions.We used about 300 000 CHAMP retrievals covering high, medium and low solar activity periods from April 2001 to August 2008.About 60 000 GRACE retrievals were used within the time period of April 2008 to December 2010.The IRO data from CHAMP and GRACE are processed by the German Aerospace Center (DLR) Neustrelitz and available at http: //swaciweb.dlr.de/ for registered users.Additionally, we used about 2.5 millions hmF2 estimates reconstructed from COS-MIC IRO observation.The COSMIC IRO data are routinely processed by the COSMIC Data Analysis and Archive Center (CDAAC) and available at http://cosmic-io.cosmic.ucar.edu/cdaac/index.html.
The IRO technique allows all-day all-season monitoring of the Earth's atmosphere (Hajj and Romans, 1998;Jakowski et al., 2002).Typical polar orbits of the LEO satellites together with the daily rotation of the Earth extend the data coverage over the globe.Thus, IRO data include day and night, summer, winter and equinoxes at high, medium and low latitudes.In general, the retrieved NmF2 and hmF2 by IRO techniques are in good agreement with ionosonde observations at medium latitudes (Jakowski et al., 2004;Angling, 2008).However, the accuracy of the reconstructed electron density profiles may degrade as a consequence of strong spatial gradients in the ionosphere, especially in the equatorial anomaly regions (Yue et al., 2010).
A large number of ionosonde data from three different sources, namely Space Physics Interactive Data Resource (SPIDR), Ionosphere Prediction Service (IPS), and National Oceanic and Atmospheric Administration (NOAA) was used in the present study.The National Geophysical Data Center (NGDC) of the United States of America (USA) provides historical and present ionosonde data records such as foF2, hmF2, foE etc. to the scientific community via SPIDR.
The SPIDR database (available at http://spidr.ngdc.noaa.gov/spidr/) currently contains over 60 years of ionospheric data from over 200 ionosondes worldwide.However, we found that for many stations hmF2 data are completely missing and for some stations data are available only for short periods.
Using the available data we computed long-time median and mean of the solar cycle variation, annual-semiannual variation and diurnal variation of hmF2 for individual stations.In some cases the median and mean variations do not follow the same pattern or largely deviate from the typical diurnal pattern, such as high hmF2 values during nighttime, sharp decrease in morning hours and gradual increase during daytime.In such cases the station is excluded from the database.Finally, we are able to use data from only 37 ionosonde stations from SPIDR.We sorted hmF2 data for selected stations and used medians for further processing of the data.
The IPS Australia is another source of historical ionospheric data used in this study (available at ftp://ftp.ips.gov.au/wdc-data/iondata/au/). The IPS stations are welldistributed over Australia.The IPS does not provide hmF2.Instead, it provides the propagation factor M(3000)F2.We used M(3000)F2 from 27 ionosonde stations provided by the IPS.Different techniques for estimating hmF2 from M(3000)F2 are reviewed in Dudeney (1983).However, we follow the algorithm/approach used in the NeQuick model, which was originally given by Radicella and Zhang (1995) based on the Dudeney (1978Dudeney ( , 1983) ) formula.For readers of this paper, we include the algorithm in the following: in which In Eq. ( 1), hmF2 is in km where M = M(3000)F2 and foF2 is the critical ionosonde frequency related to the peak electron density by NmF2 = 1.24×10 −2 (foF2) 2 .The quantity foE is the critical frequency of the ionospheric E layer.Like the NeQuick model, for computing foE we follow the Titheridge model (Leitinger et al., 1995;Titheridge, 1996), which is based on the seasonal relationship of foE with the solar radio flux F10.7 and solar zenith angle χ.
where a e is a seasonal term and its values are given in Nava et al. (2008) and χ eff is the effective solar zenith angle defined by Using Eqs. ( 1)-( 6) we converted M(3000)F2 to corresponding hmF2 values.Then, we computed the solar cycle variation, annual-semiannual variation and daily variation of hmF2 for individual stations.Depending on the mean and median analysis of hmF2 variation, we excluded stations for which the median and mean variations largely deviate from each other.Additionally, we used M(3000)F2 data from five ionosonde stations provided by NOAA via archives at ftp:// ftp.ngdc.noaa.gov/ionosonde/data/.As before, we converted M(3000)F2 to corresponding hmF2 values using Eqs.( 1)-( 6) and checked them before using in further processing.To view the ionosonde data coverage, a global map of used ionosonde stations from SPIDR, IPS and NOAA sources is given in Fig. 1.Fifteen verification stations that were used for validation purposes are also indicated in the map.We arranged a large hmF2 database, bringing together IRO data and ionosonde data.The database includes different combinations of data, like day and night, summer, winter and equinoxes, high, medium and low solar and geomagnetic activity conditions, and high, medium and low geographic/geomagnetic latitudes.After completion, we sorted and thus minimized the database to reduce the computational complexity in the fitting procedures.We sorted the database with respect to F10.7 variation, seasonal variation (i.e.day of year), local time variation and geomagnetic latitude and longitude variations.To consider the seasonal variation, hmF2 values were averaged for 27 day-intervals and the 14th day was taken as the reference day.The spatial resolution in the meridional direction was limited to 2.5 • .In the zonal direction, the maximum spatial resolution was 5 • at the equator and the resolution was gradually decreased to 360 • at the poles.The local time resolution was limited to 1 h.Thus, the length of the input data matrix was reduced to about 1 million values.
where the factors F 1−4 contain explicitly the model functions including model coefficients that describe four main dependencies of hmF2.The factor F 1 describes the daily variation with local time (LT in hours) as where V D , V SD and V TD are the angular phases of the diurnal, semi-diurnal and terdiurnal harmonic components, respectively.The functions cosχ * and cosχ * * describe the dependency on the solar zenith angle χ as where φ is the geographic latitude and δ is the declination of the sun.The value P F 1 = 0.4 in Eq. ( 11) is chosen in such a way that the term cosχ * * has always a positive contribution.
The factor F 2 describes the annual (A) and the semi-annual (SA) variation of hmF2 as given by Eq. ( 12): in which The phase shifts are best fitted as doy A = 181 days and doy SA = 49 days for the annual and the semi-annual variation, respectively.There is a geomagnetic control over the structure of the peak electron density and its height.Therefore, the peak height depends on the geomagnetic latitude φ m .For simplicity, we used a simple dipole representation of the Earth's magnetic field instead of using any multi-pole representation, such as the International Geomagnetic Reference Field (IGRF) model (Mandea and Macmillan, 2000).The latitudinal distribution of hmF2 shows a maximum at the geomagnetic equator, gradually decreasing on both sides of the equator.Our investigation shows that the peak over the geomagnetic equator is prominent during daytime, but becomes weaker during nighttime.Considering this, we modelled latitudinal distribution of hmF2 in connection with the local time variation in such a way that the magnitude of hmF2 peak is maximum at 14:00 LT and minimum during nighttime.Thus, the latitudinal distribution of hmF2 is modelled by the following expression: The half widths of the Gaussian function are best fitted as σ φ1 = 40 • , σ φ2 = 20 • and σ LT = 4 h.We have found a strong solar activity dependence of hmF2, which is modelled by the following equation: where F10.7 is the solar radio flux commonly measured in flux units (1flux unit 10 −22 W m −2 Hz −1 ) and δ F10.7 = 10.8 flux units.Equations ( 7)-( 15) explicitly describe the functional dependencies of hmF2 on local time, season, geographic and geomagnetic latitudes, and solar cycle variations.The equations contain 13 unknown polynomial coefficients in addition to a few empirically fixed known parameters.In the next section, we derive the polynomial coefficients by applying the model approach to the model database.

Modelling results
In the previous section, we formulated a set of nonlinear equations that explicitly contain model functions and coefficients.Now, we consider that the model coefficients and observation data are related through a nonlinear system of equations.Then, by non-linear least square methods, we obtained a set of 13 model coefficients.The coefficients best fit the data in the sense of minimizing the sum of squares of residual errors.To assess the degree of certainty for model coefficients, we computed the approximate covariance matrix.The standard deviations of the individual model coefficients were estimated by taking square roots of the diagonal elements of the covariance matrix.The solution coefficients and their percentage standard deviations are given in Table 1.The estimated standard deviations (STD) confirm that the model coefficients have a large degree of certainty.
It should be noted that the coefficients are related to the input data set and they may change if another or additional extended data set is used.The same is true for all the parameters fixed at certain values in the previous section, i.e. they may change when other data sets are used.Additionally, to examine the reliability of the IRO data for ionospheric parameters modelling, we computed the 13 model coefficients using only IRO data in the database.The computed model coefficients are also given in the Table 1.However, the parameters fixed at certain values in the previous section are kept the same.It should be mentioned that all comparisons and validations of NPHM given in this paper are done using coefficients obtained for IRO + ionosonde data.The only exception is Fig. 9 where both sets (IRO + ionosonde and IRO only) of coefficients are used to compare their relative performance and verify the reliability of the IRO data in ionosphere modelling.
To assess the overall fitting quality, we computed model residuals, i.e. the differences between the input data and model values.Then, we computed percentage residuals by (100 × model residual/input) %.The histogram of percentage residuals is plotted in Fig. 2. The mean deviation, standard deviation (STD) and root mean squared (RMS) estimates of percentage residuals are calculated (see in Fig. 2).
We see that the histogram is normally distributed with a small bias of about −1 % and shows no obvious asymmetric pattern.The RMS and standard deviation are found to be equal with a value of about 10 % each.Some example plots of model results and input data as a function of local time, day of year (doy) and geomagnetic latitude are given in Figs.3-5 Figure 3 displays sample plots for the local time variation of hmF2 during a summer and a winter day at high, mid and low latitudes at 0 • meridian and F10.7 = 80 flux units.We see that the model follows the data trend for most cases.In Fig. 4, the seasonal variations during daytime and nighttime are plotted at mid-latitudes for the same solar activity level.We see that during nighttime the model follows the data trend, whereas during daytime it overestimates hmF2 compared to the input data.Figure 5 gives sample plots for the geomagnetic latitude dependency during daytime and nighttime.We see that during daytime the model can successfully follow the input data.During nighttime, it tries to follow the data trend although the data trend is not clear in the input data.

Validation
For validation purposes, we did comprehensive test comparisons with observational IRO data and ionosonde data from selected time periods, and also with hmF2 estimations from the NeQuick model.It should be noted that the data used for validation purposes are already used as input data sets in the adjustment procedures.However, using a large input database (about 1 million input data sets) we derived only 13 model coefficients; therefore, we assume that it may not be unjustified to use some input data for validation purposes as well.
For comparisons to observational IRO data, we used a large database that includes five years (2002)(2003)(2004)(2005)(2006) of CHAMP data and three years (2007)(2008)(2009) of COSMIC data.The IRO data include data from day and night at high, medium and low latitudes.The data contain a mix of high and low geomagnetic activity conditions.Since we are concerned about the ionospheric F2 layer peak height, the observations that exceed the limit 200 km < hmF2 < 550 km are excluded from comparisons.The hmF2 is calculated by NPHM as well as by NeQuick model at the same location and time window as the IRO data.The percentage residuals are then computed by 100 × (obs − model)/obs % for each year of data and their histograms are shown in Fig. 6.
The corresponding RMS, mean and standard deviation of percentage residuals are given in Fig. 6, and also summa-rized in Table 2.The table shows that both the NPHM and NeQuick residuals have negative biases for all years considered.Comparing NPHM and NeQuick models, we see that the biases are less for NPHM for all cases.The maximum biases (absolute values) are found −4.8 % and −10.model for all years.The maximum STDs are found 13 % and 16.2 % for NPHM and NeQuick models, respectively, during 2002 and the minimum STDs are found 10.9 % and 13.7 %, respectively, during 2009.In general for both models, the RMS and STD of residuals are slightly higher during times of high solar activity compared to those during low solar activity periods.Comparing the statistical estimates given in Table 2, we see that although the NPHM performs slightly better than the NeQuick model, their differences are not much.
For validation with ionosonde data, we chose two onemonth periods from high and low solar activity in years 2002 and 2006, respectively.We chose the months May for 2002 and December for 2006; they correspond to the Northern Hemisphere summer and winter, respectively.Depending on the data availability during the specified one-month period, we selected 15 reference ionosonde stations distributed on both sides of the equator at high, medium and low latitudes covering the American, European and Australian longitude sectors.The solar radio flux index F10.7 is the main driving function of the proposed NPHM.However, the NPHM is climatological, i.e. it maps the long-time average behaviour of the peak density height.The small-scale features of the peak density height are smoothed out in the averaging and fitting procedures.It should be mentioned that an empirical model based on climatology cannot predict the actual variability and dynamics of the ionosphere.However, in operational ionospheric parameters reconstruction using real-time observations, the use of any empirical model as a background model is very helpful (see use of the background model for TEC reconstruction, Jakowski et al., 2011b).
Our comprehensive validation studies using CHAMP and COSMIC IRO observations and ionosonde measurements show that the performance of the new model is very similar to that of the NeQuick model.The NeQuick as well as the IRI use the CCIR M(3000)F2 maps for hmF2 estimations at any place and time.The IRI uses modified Bilitza et al. (1979) equations for M(3000)F2 conversion to hmF2 whereas the NeQuick model uses the algorithm given by Radicella and Zhang (1995) based on the Dudeney (1978Dudeney ( , 1983) ) formula.In the present paper we did not compare our model with the IRI.However, since both the NeQuick and IRI use CCIR maps for hmF2 estimations, it is expected that comparisons with the IRI will be quite similar with those with the NeQuick model.
For global distribution of hmF2, our new model uses only 13 polynomial coefficients.Thus the number of coefficients measurements.Comparisons between NPHM and electron density NeQuick models for RMS estimates of their differences from observational data show that during high solar activity period the RMS deviations are about 13 % and 18 % for NPHM and NeQuick models, respectively.During low solar activity periods, the corresponding RMS estimates are 12 % and 16 %, respectively.The performance of the new model may be further improved in the future by extending the model database by integrating available topside sounder data and utilizing IRO data from upcoming satellite missions.
Due to lack of observational data for the specified onemonth period May 2002, we selected May 2001 for Rome and May 2003 for Puerto Rico as the test periods.The station name, location and test period are given in Table 3. Figures 7 and 8 compare the NPHM and NeQuick model results with the ionosonde measurements as a function of Universal Time (UT) at the selected station locations.The hmF2 values are averaged at each UT hour for all days in May and December.Figures 7 and 8 show the comparisons to the observational ionosonde data for the Northern and Southern Hemisphere stations, respectively.

Figure 7
shows that NPHM performs better than the NeQuick model at Tromsø, Juliusruh, Rome, Dyess and Puerto Rico during the selected high solar activity month May.Only at Eglin AFB NeQuick performs better than NPHM.During low solar activity period December 2006, NPHM performs better at Dyess, Eglin AFB and Kwajalein, whereas NeQuick performs better at Tromsø, Juliusruh and Athens.Figure 8 shows that at Southern Hemisphere stations Hobart and Macquarie Island, NPHM performs better

Figure 7 .Fig. 7 .Figure 8 . 4 Fig. 8 .
Figure 7. Monthly mean of NPHM and NeQuick as a function of UT in comparisons with 2 ionosonde observation at northern hemisphere stations.3 Fig. 7. Monthly mean of NPHM and NeQuick as a function of UT in comparison with ionosonde observation at Northern Hemisphere stations.

www.ann-geophys.net/30/797/2012/ Ann. Geophys., 30, 797-809, 2012 800 M. M. Hoque and N. Jakowski: A new global model for the ionospheric F2 peak height 3 Modelling approach
In our recent papers(JakowsHoque and Jakowski, 2011b)and Jakowski, 2011b), we developed basic approaches for modelling ionospheric parameters on global scale, e.g.NTCM (Neustrelitz TeC Model) and NPDM (Neustrelitz Peak Density Model) approaches.Following the same basic approach, in the present work we developed a set of nonlinear equations describing the dependencies of hmF2 on local time, season, geomagnetic field and solar activity as

Table 1 .
Optimal set of model coefficients and their percentage standard deviations.

Table 2 .
Estimates of percentage residuals for NPHM and NeQuick model in comparison with IRO data.
than the NeQuick during May 2002.The NeQuick performs better at Vanimo, Townsville and Port Stanley.During December 2006, at Hobart and Mawson NPHM performs better, whereas at Vanimo, Jicamarca and Townswille NeQuick model shows comparatively better performance.Figure 9 also confirms these findings.As already mentioned, we determined two sets of model coefficients (see Table1) depending on the used database: