© A.V. Mikhailov et al., Published by EDP Sciences 2014

1. Introduction

There is a specific need to develop methods for monitoring the thermosphere, in order to allow satellite systems to compensate for atmospheric drag effects if a booster is available. Global 3D first-principle models like the Thermosphere Ionosphere Electrodynamics General Circulation Model (Richmond et al. 1992) or the Coupled Magnetosphere Ionosphere Thermosphere (CMIT) model (Wang et al. 2004) are powerful research tools in our efforts to understand the physical processes in the Earth’s environment. However due to the extremely complicated mechanisms of the upper atmosphere interaction with layers from above and below, problems with the models’ performance under disturbed conditions have been stressed repeatedly. Many of these problems have been attributed to the lack of reliable estimations of the ionospheric drivers. Another group of broadly used models is the semi-empirical atmospheric density models, such as DTM-94 (Berger et al. 1998), DTM-2000 (Bruinsma & Thuillier 2003), and MSIS family of models (e.g., Hedin 1987). Lathuillère et al. (2008) presented a comparison between CHAMP/STAR observations with the MSISE-00 model predictions (Picone et al. 2002) for disturbed conditions. This comparison has shown that “the model statistically underestimates the density disturbance by about 50% for magnetic activity levels less than Kp ~ 6. During more disturbed conditions the model can underestimate the density disturbance amplitude by more than 200%”. But such periods are the most important and interesting ones from the practical point of view keeping in mind that during disturbed conditions the atmosphere heats and expands, exerting an increased drag on satellites orbiting around the Earth. The increased atmospheric drag causes satellites to loose altitude and finally enter the dense atmosphere. Therefore the development of methods to monitor and reliably forecast atmospheric drag effects on satellites under all geophysical conditions is a very actual problem with strong economic consequences.

Mikhailov et al. (2012) proposed a new approach to monitor the thermosphere (referred to as M2012 thereafter), based on quite different principles. The Earth’s ionosphere is embedded in the neutral thermosphere (even in the F2-layer maximum the ion-to-neutrals ratio is ~ 10⁻³) and it closely follows the state of the thermosphere at least during sunlit hours. Therefore the ionosphere may serve as an indicator of the state of the thermosphere under various geophysical conditions. For this reason, the idea to extract thermospheric parameters from ionospheric observations seems to be promising. To test its feasibility at middle latitudes, the method was implemented using as input Millstone Hill ISR electron density profiles and the results were compared with CHAMP neutral gas density observations both under solar maximum (October 2002) and minimum (January 2008) activity. It was shown that the retrieved neutral gas densities were in the ±(10–15)% corridor that corresponds to the announced absolute inaccuracy of CHAMP observations. In an attempt to assess the operational potential of the method, we have used as input Digisonde automatically scaled EDP, extrapolated with the NeQuick profiler. As it was expected, because of the ionograms autoscaling errors reflected in the input EDP and the assumptions related to the topside extrapolation, the accuracy of the results is degraded comparing to what we got using ISR EDP, but nonetheless, there was a general agreement between CHAMP and the extracted neutral densities from M2012, especially under solar maximum conditions. It was stressed that mid-latitude routine Digisonde Ne(h) can be considered as input parameters in the proposed method for routine retrieval of thermospheric parameters after careful reprocessing to minimize autoscaling errors. Similar conclusions were obtained for the equatorial F2-region (Mikhailov et al. 2013) where ISR and Digisonde EDP from Jicamarca Observatory were used for the retrieval of thermospheric parameters. In this case, CHAMP and GRACE neutral gas density observations were used to evaluate the results of the method.

In this paper we examine the performance of the M2012 method, using as input an alternative source of EDP. In particular, we use Ionospheric Radio Occultation (IRO) measurements on board LEO (Low Earth Orbit) satellites. Radio Occultation (RO) is one of the few techniques that can probe both top and bottom sides of the ionosphere. The technique for sounding the Earth’s atmosphere was demonstrated by the proof-of-concept GPS Meteorology (GPS/MET) experiment in 1995–1997 (Ware et al. 1996). Following GPS/MET, additional missions, that is, the Challenging Minisatellite Payload (CHAMP; Wickert et al. 2001) and the Satellite de Aplicaciones Cientificas – C (SAC-C; Hajj et al. 2004), have confirmed the potential of RO sounding of the ionosphere, stratosphere, and troposphere. The FORMOSAT-3/COSMIC (Rocken et al. 2000) also applies this technique, with primary scientific goal to demonstrate the value of near-real-time RO observations in operational numerical weather prediction, including space weather science and operations. Having access to a wealth of IRO profiles, we assess the performance of the thermospheric retrieval method M2012 using as input IRO EDP obtained from CHAMP and FORMOSAT-3/COSMIC missions.

Unlike Digisonde or ISR observations which provide local overhead Ne(h) distributions, IRO observations have a number of constrains, due to their retrieval method, that have to be considered when assessing the results obtained from the M2012 method.

In summary, the aims of the paper are formulated as follows:

1. Testing the method’s performance using simultaneous IRO Ne(h) observations and CHAMP neutral gas density measurements.

2. Assessing the results taking into account the limitations of the IRO observing technique.

3. Evaluating whether the developed method can be used in practice to retrieve thermospheric parameters (neutral composition and temperature) from routine IRO EDPs.

In the following Section 2 we present the data used for this analysis, the observing method, the limitations and accuracy. The method’s results with IRO profiles are provided in Section 3 during solar maximum and solar minimum conditions separately and a comparison with CHAMP neutral density profiles is attempted. In Section 4 we evaluate the results in respect to the accuracy and limitations of IRO EDPs. Finally in Section 5 we discuss the results and summarize the main conclusions.

2. Data presentation

2.1. Ionospheric radio occultation measurements

GPS RO is a space-borne remote sensing technique providing accurate, all-weather, high vertical resolution profiles of atmospheric parameters (Melbourne et al. 1994; Ware et al. 1996; Kursinski et al. 1997; Rocken et al. 1997; Wickert et al. 2001; Hajj et al. 2002).

For the purpose of this analysis we use IRO data from CHAMP and COSMIC satellites.

CHAMP satellite was launched into a near polar orbit (I = 87°, h = 450 km) by a Russian COSMOS rocket on 15 July 2000 (Reigber et al. 2000). Due to the high inclination, the radio occultation measurements obtained onboard CHAMP cover the global ionosphere. However only IRO measurements from middle latitudes have been used in this work because the assumption that the ionospheric flow follows the thermospheric flow is not generally valid at high latitudes under disturbed conditions. Moreover it is expected that the middle latitude IRO profiles are more reliable in a comparison to high-latitude ones, due to strong spatial irregularities in the high latitude ionosphere. The satellite is equipped with a dual frequency GPS receiver that enables the analysis of the 0.1 Hz sampled navigation data as well as GPS radio occultation measurements. Thus, the actual state of the ionosphere is permanently monitored near the CHAMP orbit plane.

The advanced “Black Jack” GPS receiver developed by the Jet Propulsion Laboratory (JPL) provides, on one hand, precise time and orbit information. On the other hand, in the radio occultation or limb sounding mode, the receiver measures the GPS carrier phases starting at CHAMP orbit tangential heights down to the Earth surface with a sampling rate of 1 Hz. Furthermore, the 0.1 Hz sampling navigation data may be used for reconstructing the electron density distribution in the topside ionosphere by data assimilation (Heise et al. 2002).

From about 200 IRO measurements about 150 EDPs are successfully retrieved per day. Because the processing system works automatically, EDP outliers cannot be avoided, but the number of such “unrealistic” profile outliers is in principle less than 1% (Jakowski et al. 2004).

FORMOSAT-3/COSMIC mission (C/F3) consists of six microsatellites that were launched on 15 April 2006 and reached their mission orbit of 800 km in December 2007, with a separation angle between neighboring orbital planes of 30° longitude. Their orbital configuration gives global coverage of approximately 2000 soundings per day, distributed nearly uniformly in local solar time. Today, more than 2,700,000 EDP are available at http://cdaac-ftp.cosmic.ucar.edu/cdaac/ for various geophysical conditions.

The C/F3 EDP is retrieved from GPS-RO data along the GPS-LEO radio links near the ray path tangent points. Recent studies show that ionospheric electron densities derived from the RO sounding around and above the F2 peak are reasonably accurate, whereas those below the F region should be used with great caution because of the assumption of spherical symmetry used in the Abel inversion (Lei et al. 2007; Kelley et al. 2009; Liu et al. 2010; Yue et al. 2010). As a result, Lee et al. (2012) apply a C/F3 EDP observation error assumed to be the sum of a 10% instrumentation error and the estimated Abel inversion error percentage (Liu et al. 2010; Yue et al. 2010).

2.2. CHAMP neutral gas density observations

Retrieved neutral gas density was compared to CHAMP/STAR observations. Detailed information on these observations as well as the accuracy estimates can be found in Bruinsma et al. (2004, 2006) as well as in the Sutton (2008) thesis. For our comparison two issues are important: (i) what is the uncertainty in the observed absolute density and (ii) how much is the contribution of He to the total neutral density at the heights of the comparison. The latter is important as we retrieve only [O], [O₂], and [N₂] concentrations, but the observed neutral density includes [He] as well. For CHAMP observations “the total uncertainty in the absolute density value is at the 10–15% level for up to moderate geomagnetic activity conditions (Ap = 15)” (Bruinsma et al. 2004). The majority of our analyzed cases are under this condition.

The CHAMP neutral gas density observations used in this paper were downloaded from http://sisko.colorado.edu/sutton/data.html. The observed neutral densities in the vicinity of the RO tangent points (the difference in longitude and latitude ≤5°) were reduced to the location and UT of RO tangent points using the MSISE-00 (Picone et al. 2002) thermospheric model and the following expression:where ρ _CHAMP is the total neutral gas density measured with CHAMP, MSIS_tp is the total neutral gas density obtained from MSIS model in the tangent point tp, MSIS_CHAMP is the total neutral gas density obtained from MSIS model in the CHAMP position.

The height for the reduction is selected to be close (the difference <5 km) to the height of CHAMP neutral gas density observations to minimize possible errors due to the MSISE-00 imperfectness.

The retrieved neutral gas densities should be corrected for [He] to be compared with the observed and modeled values. The contribution of He to the total density is not essential under high solar activity and CHAMP observations (heights ~400 km), but this contribution may be noticeable under solar minimum in 2007–2008. We could confine ourselves with the MSISE-00 empirical model, which provides all neutral species concentrations as the output parameters, but some model comparisons with CHAMP and GRACE neutral densities have shown the advantages of the JB-2008 model over the MSISE-00 one (Bowman et al. 2008, their Fig. 10; Shim et al. 2012). Therefore, it would be interesting to estimate the [He] contribution to the total neutral density with the JB-2008 model as well. The JB-2008 model is known to give only neutral temperature and density height profiles. Therefore the modeled neutral gas density should be split to individual neutral species [O], [O₂], [N₂], and [He]. This can be done supposing a barometric height distribution for individual species (Mikhailov et al. 2013). The contribution of [He] to the total modeled (MSISE-00 and JB-2008) neutral gas density turned out to be less than 2% at the height of a comparison for solar minimum in 2008. Similar conclusion was obtained earlier (Mikhailov et al. 2013). This contribution is much less than the uncertainty in ρ observations (10–15%) we are working with. Although the retrieved ρ values were corrected for the He contribution using the model ratio ρ _tot/ρ _tot-He at the height of a comparison, this does not appreciably affect our comparison with the observations.

Unlike ionosonde or ISR observations that provide us with Ne(h) variations over time at the same geographic location, in the case of IRO observations we have one “stand-alone” Ne(h) profile without any prehistory, therefore only a stationary version of the M2012 can be applied.

3. Testing of the method’s performance

The method was tested both for solar minimum (2007–2008) and solar maximum (2002) conditions. Only CHAMP IRO and neutral gas density observations are available for solar maximum while COSMIC IRO and CHAMP neutral gas density measurements were used under solar minimum conditions.

3.1. Solar minimum conditions

The method was tested for 87 cases under solar minimum conditions and the results are provided in tabular form in the Appendix. For 37 cases no solution could be obtained (Table A3). In the following we focus on the 50 cases for which solutions could be obtained. To classify these cases, we calculated the distribution of the R = ρ _cal/ρ _obs ratio and then we estimated the average R _ave and the standard deviation SD. Here ρ _cal is the neutral density extracted with the M2012 method and ρ _obs are the corresponding CHAMP observations. Figure 1 gives the histogram of R. The dashed vertical lines correspond to R _ave ± SD. The majority of cases (77%) fall within the range of R _ave – SD ≤ R ≤ R _ave + SD. These cases are listed in Table A1, while those residing outside this range (33%) are listed in Table A2. Although IRO EDP are from both hemispheres in the sample of cases with R _ave – SD ≤ R ≤ R _ave + SD, the majority of the cases considered are found in summer and equinoctial seasons. Geomagnetic activity was either low or moderately elevated.

Fig. 1. Distributions of R = ρ cal/ρ obs ratio for the retrieved cases with the M2012 method (top panel), and the corresponding distribution based on JB-2008 and MSISE00 neutral gas densities models for solar minimum conditions (middle and bottom panel). The average R ave and the standard deviation SD values are given. The area between dashed lines corresponds to R ave ± SD.

The retrieved neutral gas densities are compared to JB-2008 and MSISE00 model values. Comparing the R _ave obtained for the three distributions presented in Figure 1, the minimum R _ave value (average ρ _cal/ρ _obs is closer to 1.0) takes place with the M2012 method, although a long tail of low ρ _cal/ρ _obs is seen.

Our previous analysis (M2012) has shown that low retrieved neutral gas densities are related to low observed hmF2 values. The underestimation of hmF2 in IRO observations will be discussed in the next section.

Focusing on the cases within the R _ave ± SD range, we provide some statistical metrics (i.e., mean relative deviation MRD and the bias with respect to the observed values) for a comparison between the M2012 and the performance of empirical models. The M2012 method gives MRD = 9.3% and the bias = −0.07 × 10⁻¹⁵ g cm⁻³; JB-2008 model gives MRD = 16.4% and the bias = 0.39 × 10⁻¹⁵ g cm⁻³; MSISE-00 model gives MRD = 28.9% and the bias = 0.82 × 10⁻¹⁵ g cm⁻³. This test demonstrates that the proposed method tends to have better accuracy than modern empirical models in those cases when an acceptable solution can be obtained with the M2012 method. It should be stressed that the retrieved neutral gas densities are close to the announced absolute uncertainty ±(10–15%) of the neutral gas density observations with the CHAMP satellite (Bruinsma et al. 2004). For a quick visual inspection, the plot of the modeled vs. observed neutral gas densities extracted from the cases that fall within the R _ave ± SD range is given in Figure 2. This graphical representation and the statistical results show that the retrieved gas densities are more centered with respect to the observations compared to modeled values which are essentially biased overestimating the observations. However, as it was mentioned earlier not all IRO profiles can provide an acceptable solution. There is a class of EDP which formally gives solutions but the retrieved neutral gas densities turn out to be far from the observed values. These are the cases outside the R _ave ± SD range of the distribution (Fig. 1). Overall 11 such cases encountered during our analysis and they are listed in Table A2.

Fig. 2. Retrieved and modeled neutral gas densities from Table A1 are given vs. the observed values. Note that the retrieved densities are more centered with respect to the observed ones while modeled values are biased overestimating the observations. MRD and bias are given for a comparison.

Considering all cases for which a solution could be obtained we have recalculated the statistical metrics. In this case the M2012 method gives MRD = 15.9% and the bias = −0.20 × 10⁻¹⁵ g cm⁻³; JB-2008 model gives MRD = 17.0% and the bias = 0.41 × 10⁻¹⁵ g cm⁻³; MSISE-00 model gives MRD = 28.0% and the bias = 0.80 × 10⁻¹⁵ g cm⁻³. It is seen that the MRD of our method is the least and is still close to the announced absolute inaccuracy of CHAMP neutral gas density observations. The least bias tells us that the retrieved values are more centered with respect to the observations compared to the empirical models. This result is important from practical point of view: if a solution exists then we may expect that the retrieved neutral gas densities on average should be close to the real values, the results being better than the empirical models provide. Table A2 shows that, for the retrieved densities large deviations may take place, but large deviations may be obtained in respect to CHAMP observed densities, for the empirical models as well.

Concerning the third class of IRO EDP for which no solutions could be found, overall 37 such cases were encountered during our analysis of solar minimum conditions, listed in Table A3. Analysis of the numerical values shows that these cases are not related either to solar or to geomagnetic activity. The feature that tends to characterize these cases is the unrealistic hmF2 values. The majority of these cases manifest very low hmF2, although cases with large hmF2 are also among them (e.g., see 01.01.07; 31.01.08; 12.07.08 in Table A3). One can also find cases with unrealistic large NmF2 (e.g., see 01.01.07; 29.01.08; 29.10.08 in Table A3). In Section 4 “Assessment of the results” it will be shown that for these cases an inconsistency is observed between hmF2 and NmF2 and for this reason such IRO EDP are not consistent within the mid-latitude daytime F2-layer formalism adopted by the M2012 method.

Summarizing, it is concluded that about half of the analyzed IRO EDP under solar minimum could be developed with the M2012 method to result in neutral gas density close to CHAMP observations. The results were shown to be, in general of higher accuracy than those obtained with the JB-2008 and MSISE-00 empirical models, comparing to the corresponding neutral density observed values. However, more than half of the analyzed IRO profiles could not be developed at all. Insufficient accuracy of such IRO EDP is considered as a possible reason. This problem will be analyzed in the following section. Finally, it should be noted that within the same day we may have IRO EDP profiles that give a solution within the R _ave ± SD range, or outside this range, or cases that do not provide a solution at all.

3.2. Solar maximum conditions

Now let us consider testing results for solar maximum conditions using the observations for 2002. Only CHAMP both IRO EDP and neutral gas density observations are available for solar maximum and the total number of cases analyzed is limited to 36. An additional difficulty related to CHAMP IRO observations is that the electron density is available only below 400–410 km and the upper boundary condition is taken at 400 km height. This is too close to the F2-layer peak which may occur at ~(320–330) km and this may affect the fitting process resulting in incorrect solutions.

Both hemispheres were considered, but again summer and equinoctial cases dominate. Unlike solar minimum conditions, the list now includes dates with strong geomagnetic disturbances up to Ap = 70. During solar maximum conditions, a solution could be obtained in all 36 cases considered.

By analogy with the analysis of cases under solar minimum we have drawn in Figure 3 (top panel) the histogram for the R = ρ _cal/ρ _obs ratio, then found an average R _ave = 1.20 and SD = 0.29. The corresponding results obtained using JB-2008 and MSISE00 models are plotted in histogram format as well, in the middle and bottom panels of Figure 3.

Fig. 3. Same as Figure 1 but for solar maximum conditions.

Considering only the cases inside the R _ave ± SD range, we calculated the mean relative deviations (MRD) and bias to quantify the comparison with observations. The proposed method gives a bias of 0.72 g cm⁻³ and MRD = 15.5%; JB-2008 model gives a bias of 0.27 g cm⁻³ and MRD = 12.0%; MSISE-00 model gives a bias of −0.24 g cm⁻³ and MRD = 9.3%. Unlike the results for solar minimum, the retrieved gas densities manifest worse accuracy than the empirical models provide. This contradicts the results obtained for Millstone Hill ISR EDP with the M2012 method where the efficiency of the proposed method was proved to be higher for solar maximum activity compared to solar minimum. The retrieved densities demonstrate a positive bias which is seen from Figure 3 and Table A4 (the last column) – the retrieved ρ are systematically larger than the observed values. These deficiencies will be further discussed in the Discussion part.

Similarly to solar minimum conditions, there are cases with deviations outside the R _ave – SD ≤ R ≤ R _ave + SD range listed in a tabular form in the Appendix (Table A5). Again these cases are not related to enhanced geomagnetic activity and take place both under quiet and disturbed conditions and may occur for the same dates when we have cases that provide a solution within the R _ave ± SD range. In general all the solutions found for solar maximum are not successful enough. This may be related to the accuracy of the IRO EDP as it will be discussed later.

4. Assessment of the results

In order to assess the results obtained from our method, and especially analyze the cases that did not result in a solution, we will analyze how the precision of IRO EDP affect the performance of the M2012 method. The precision of IRO EDP for both CHAMP and COSMIC observations has been determined in various studies in the past mainly through comparisons between IRO EDP with ISR and ionosonde data.

Early studies by Hajj & Romans (1998) presented a comparison between IRO EDP from the GPS/MET experiment with one Millstone Hill ISR Ne(h) profile, as well as with NmF2 observations obtained from a global network of ionosondes during a 20-day period in April–July 1995. Their RO/ISR comparison looks fairly good especially for ISR observations with 320-μs pulse mode in the whole 200–500 km height range. The RO NmF2 observations agreed with the ionosonde measurements with a deviation within 20% (10% in foF2) and they were proven to be essentially unbiased. However, it should be noted that GPS occultation took place within ≈1100 km radius (10 degrees) around the ionosonde stations and this is larger than the correlation length for solar minimum conditions (Kiseleva et al. 1971; McNamara 2009). Therefore, formally not all observations could be used for a comparison. Schreiner et al. (1999) found a 13% (~0.5 MHz) deviation of GPS/MET measurements with respect to foF2 ionosonde observations (45 ionosondes, February 20–23, 1997). One of the first comparisons between CHAMP RO measurements with in situ electron density PLP observations as well as with ionosonde foF2 and hmF2 data was performed by Jakowski et al. (2002). Langmuir probe data obtained on board the CHAMP satellite were directly compared to RO observations during 26 days in 2001. The horizontal distance between Planar Langmuir Probe (PLP) and RO measurements was less than 700 km. The results demonstrate a quite consistent correlation between the two datasets. On the other hand, about 100 measurement coincidences within the cross section diameter of ≈1900 km (this is larger than the spatial correlation length even for solar maximum) between RO and ionosondes were found. The bias between RO and ionosonde parameters amounted to −1.7% for foF2 and −4.1% for hmF2 indicates a slight underestimation in the RO data. The derived MRDs were found to be 17.8% for foF2 and 13.1% for hmF2. Assuming the peak height hmF2 equal to 300 km, the results correspond to ≈−12 km for the bias and ≈40 km for the average hmF2 deviation. The deviation seems to be large, but it should be kept in mind that ionosonde hmF2 were estimated using the approximate dependence by Dudeney (1983).

A more extended comparison of COSMIC RO profiles with Arecibo ISR observations in June 2006 was presented by Kelley et al. (2009). This comparison included 32 profiles in overall, which were obtained by using the Abel transform method that was developed in two versions independently at UCAR and JPL. No clear conclusion could be extracted for the nine daytime comparisons (their Fig. 7), as the agreement ranged in all possible levels. An interesting result is that the two methods of Abel transformation may give quite different EDP telling us that the data processing methods are not straightforward. Both NmF2 and hmF2 scatterplots (their Figs. 10 and 11) manifest large scatter, but NmF2 is seen to be determined better than hmF2. The SD of the former is 1 × 10⁵ cm⁻³, while no statistical results are given for hmF2. It is pointed out that the best agreement between hmF2 IRO values with the ISR ones was obtained near 300 km, but IRO values are lower than ISR values below this height and higher than ISR values above it.

Similar but more representative results can be found in Tsai et al. (2009) who made a comparison of COSMIC observations with 49 worldwide distributed ionosondes for the period of 20 June–27 September, 2006. Their Figure 6 exhibits a large scatter of COSMIC foF2 vs. ionosonde values with mean relative deviation 20% (40% in NmF2). This coincides with earlier estimates by Schreiner et al. (1999) and Tsai et al. (2001). A comparison on hmF2 has shown that “the FS3/COSMIC hmF2s do not coincide well with the ionosonde hmF2s”. Their Figure 9 shows that COSMIC hmF2s are much lower than ionosonde values.

Schreiner et al. (2007) analyzing the precision of GPS RO from FORMOSAT-3/COSMIC mission concluded: “Thus retrieved EDPs are expected to have rather poor accuracy when interpreted as actual vertical profiles”. According to Schreiner et al. (2007) the largest error in the GPS IRO retrieved EDPs is due to strong horizontal gradients disrupting the assumption of spherical symmetry. This assumption, in cases with large NmF2 values, can result in either positive or negative errors larger than 10⁵ cm⁻³ at the bottom of the retrieved profiles (Syndergaard et al. 2006). This conclusion seems to support the results by Lee et al. (2012) who have attempted to assimilate FORMOSAT-3/COSMIC GPS occultation observations into the NCAR TIEGC model. A comparison of the assimilation results with the independent ionosonde observations has shown the improvement of the primary ionospheric parameter, such as NmF2 and hmF2, however 35% of RO EDP were rejected mainly due to either outlier threshold or quality control. This high rejection rate indicates that not all RO observations can be used for physical interpretations – the conclusion that also follows from our analysis (see later in this section).

Summarizing the results of the previous comparisons between RO Ne(h) observations with both ionosonde and ISR measurements it is possible to conclude the following. The agreement between RO retrieved and ionosonde foF2 values is within 10–20% (20–40% in NmF2). This is a rather large discrepancy for our method to retrieve the thermospheric parameters from EDP. However, one should keep in mind that not all comparisons were performed for the occultations and ionosondes to be within the correlation radius for foF2. No spatial and temporal reduction was applied to the observed foF2 to compare two types of observations as this follows from the publications. The most valuable comparisons are with daytime ISR EDP, but they are very limited in number and were only made for some selected days.

Therefore we have made additional comparisons using available simultaneous Millstone Hill ISR (http://madrigal.haystack.mit.edu/madrigal/) and COSMIC RO (http://cdaac-ftp.cosmic.ucar.edu/cdaac/) observations in 2007–2008 with tangent points at the F2-layer maximum height within ≤6° in latitude and longitude of Millstone Hill ISR (42.6 N, 288.5 E). The F2-layer spatial correlation length is ≈700 km at middle latitudes under solar minimum conditions (Kiseleva et al. 1971; McNamara 2009) and this is larger than a 6° radius we used to select COSMIC observations. Tables 1 and 2 give available simultaneous Millstone Hill ISR/COSMIC NmF2 and hmF2 measurements for summer and winter months along with the results of Ne(h) comparisons in terms of the mean relative deviations estimated separately for the topside and bottom side parts. Table 3 using data from Tables 1 and 2 gives statistical results of the comparison between Millstone Hill ISR and COSMIC IRO observations. For better visualization, the scatterplots of hmF2 and NmF2 determined by Millstone Hill ISR and COSMIC RO methods are given in Figure 4. Table 3 shows that the MRD for NmF2 is 20–24% (10–12% in foF2) and this coincides with the minimal earlier mentioned estimates. However COSMIC RO NmF2 values are slightly biased by 0.58 × 10⁵ cm⁻³ with respect to ISR observations (Fig. 4).

Fig. 4. Scatterplots of hmF2 and NmF2 determined by Millstone Hill ISR vs. RO method.

Although the MRD for hmF2 is relatively small ≈10% (Table 3), this corresponds to a 22 km bias for COSMIC hmF2 observations (Fig. 4), which are systematically lower than ISR values. The same tendency was found by Kelley et al. (2009) who used the Arecibo ISR observations – COSMIC hmF2s were lower than ISR values if hmF2s were below 300 km height as this is the case in our comparisons (Fig. 4). This result also agrees with the earlier conclusion by Tsai et al. (2009) who found COSMIC hmF2s to be much lower in a comparison with the ionosonde hmF2 values.

Selected cases of COSMIC coincidences above Millstone Hill are given in Figure 5, where COSMIC IRO EDP (solid lines) are overplotted with the corresponding Millstone Hill ISR EDP (triangles). Table 3 and Figure 5 show that Ne(h) in the topside agree fairly well in the two types of measurements, but this is not so in the bottom side. Although the selected winter cases in Figure 5 manifest a fairly good coincidence in the bottom side, Table 3 shows large MRD for winter months similar to summer ones. An obvious disagreement between the two observations takes place in the Ne(h) bottom side for summer months. IRO observations give a very broad Ne(h) in the bottom side with the main ionospheric peak at low heights, as for example on July 12, 2008 at 170 km (Table 1). Partly, this may be related to a moderate magnetic disturbance (Ap = 16, Table 1). It is known that under such conditions the F1 layer is most pronounced and the formation of G-conditions (NmF1 > NmF2) is more probable. But IRO observations demonstrate the same type of Ne(h) profile for very quiet days of July 09 and July 11, 2008 (Fig. 5). Possible explanation may be related to sporadic E which is very probable during daytime summer conditions at middle latitudes. The analysis of FORMOSAT-3/COSMIC RO profiles in Rome whose latitude is close to the latitude of Millstone Hill has shown a strong influence of sporadic E on the Ne(h) distribution in the bottom side (Perrone et al. 2011). Our analysis of Millstone Hill ISR profiles for the summer dates listed in Figure 5 has shown the existence of strong underlying ionization below 150 km which should influence RO profiles. We do not work with heights below 175 km, for this reason the lower part of Ne(h) is not shown in Figure 5.

Fig. 5. Daytime Millstone Hill ISR EDP (triangles) in a comparison to COSMIC RO measurements (solid line) under solar minimum 2007–2008 conditions. Six of the best (top panels) and six of the worst (lower panels) coincidences are shown.

Although large discrepancies in the EDP bottom side are mainly observed in summer, they may take place in winter as well. Such a case is seen for example on December 11, 2007 under a moderate (Ap = 16/9, Table 1) geomagnetic disturbance. On the other hand, under a stronger disturbance with Ap = 22/18 (Table 1) on December 18, 2007 we have a fairly good agreement between the two observations (Fig. 5). One may think that this is related to the direction in which the occultation is observed.

Summarizing the results of Millstone Hill ISR/COSMIC RO Ne(h) comparison for daytime solar minimum (2007–2008) conditions (totally 35 cases) we may conclude the following: in the majority of cases the topside of EDP coincides fairly well in the two types of observations and this is important for our method as the upper boundary condition for the continuity equation is taken from the topside of the observed EDP. The coincidence in the bottom part of EDP is not good for the 40% of the analyzed cases, and this may be a serious limitation for our method. This is close to the estimate by Lee et al. (2012) who rejected about 35% of IRO EDP due to either outlier threshold or quality control.

In this paper COSMIC IRO measurements for solar minimum (2007–2008) and CHAMP IRO observations for solar maximum (2002) were used to retrieve thermospheric parameters in the middle latitude daytime F2-region. Electron density spatial gradients are supposed to be the smallest under such conditions and this should help the method. However, keeping in mind the previous analysis of the IRO EDP observations we may expect a high rejection rate in the IRO measurements. The EDP with many negative electron density values and those obtained under large (>75°) solar zenith angles were rejected at the initial stage. There are many EDP with unrealistic low or too large hmF2 for which no solution can be obtained. Some of them were obtained at low (<35°) latitudes and one may think that they bear the effects of latitudinal gradients related to the equatorial anomaly. Anyway the absence of a solution tells us that such EDP cannot be described within the classic mid-latitude daytime F2-layer formalism used in our method. The quality of the IRO profiles varies from observation to observation and within the same day and we may have profiles that provide solution with high confidence, and profiles that provide no solution at all (cf. Tables A1–A3). The observed EDP is the main input to the method and the existence of a solution depends at a large degree on its quality.

To further approach the problem of cases for which no solution can be taken we have analyzed a specific case comparing ISR EDP with IRO EDP. January 31, 2008 case is an indicative example illustrating this dependence. COSMIC IRO tangent point (46.9E; 283.9E; 18.7 UT) was not far from Millstone Hill IS radar (42.6E; 288.5E) and CHAMP neutral gas density observations were available in that area as well. The two observed EDP are shown in Figure 6. Although both observations are within the spatial correlation distance (Kiseleva et al. 1971; McNamara 2009) the two profiles look different. The IRO EDP is shifted downward while NmF2 concentrations are close in the two cases. This difference turned out to be crucial as no solution could be obtained with the COSMIC profile. On the other hand both the instantaneous ISR EDP for 18.7 UT and the average ISR EDP obtained over the (18.42–18.96) UT time interval gave solutions with the retrieved neutral gas density close to CHAMP observations ρ = 3.55 × 10⁻¹⁵ g cm⁻³ reduced to 340 km. The instantaneous ISR EDP profile gave ρ = 3.67 × 10⁻¹⁵ g cm⁻³ and the average ISR profile gave ρ = 3.45 × 10⁻¹⁵ g cm⁻³.

Fig. 6. Observed with COSMIC and Millstone Hill ISR EDP on January 31, 2008.

The downward shifting of hmF2 and the whole bottom side part of the COSMIC EDP resulted in the absence of solution with the M2012 method. Table A3 gives many of such EDP with abnormally low hmF2.

According to the theory of the F2-layer, NmF2 and hmF2 are closely related and they cannot take arbitrary values which is often the case with IRO observations, for instance, large NmF2 and very low hmF2 (Fig. 5, Table A3). Obviously this inconsistency contradicts our mid-latitude daytime F2-layer formalism and no solution can be obtained in such case.

Figure 7 gives an additional illustration. We took Millstone Hill ISR EDP and the available COSMIC mid-latitude daytime observations for the same period of July 2008. COSMIC occultation tangent points might not be close to Millstone Hill, but longitudinal variations of F2-layer parameters are not strong during daytime hours in summer at middle latitudes (Bilitza & Reinisch 2008).

Fig. 7. Observed at Millstone Hill ISR in July 2008 daytime hmF2 vs. LogNmF2 in a comparison with COSMIC cases for which the solutions are within the R ave ± SD range and cases for which no solutions were obtained.

Figure 7 gives a scatter plot of Millstone Hill hmF2 vs. LogNmF2 in a comparison with COSMIC observed values. All COSMIC hmF2-NmF2 points which did not give a solution (red dots in Figure 7) are far away from the cloud of Millstone Hill points (marked with black dots). The majority of them manifest lower hmF2 values, although some of these cases exhibit abnormally high hmF2. For some of these cases the IRO NmF2 exceeds significantly the average NmF2 observed with the ISR at Millstone Hill. The dots marked with green, correspond to the COSMIC IRO EDP which gave solutions and are seen to be closer to the Millstone Hill points, especially concerning the hmF2 values. We stress once again that ISR Ne(h) which always provide acceptable solutions should be considered as correct ones. The tendency of COSMIC EDP to demonstrate lower hmF2 and larger NmF2 in a comparison to ISR observations was shown earlier (Figs. 4 and 5).

The inconsistency of IRO profiles can be also seen in a comparison with the climatologic ionospheric IRI model. Figure 8 gives scatter plot IRO hmF2 vs. LogNmF2 for all analyzed cases in a comparison with the IRI-95 (Bilitza 1997) monthly median model values calculated for the same IRO experiments. For solar minimum similar to Figure 7 all cases which did not give a solution occupy the lower part of the plot manifesting very low hmF2 while there are no points in IRI with hmF2 ≤ 200 km. Practically all IRO points are located below 280 km, while IRI gives many points above the 280 km height. Similar results demonstrate a comparison with IRI for solar maximum condition – the majority of hmF2 IRO observations are below the IRI model values. This basic and important result has been stressed repeatedly in our paper using various comparisons. Low hmF2 as it was shown in M2012 results in low retrieved neutral gas density. This may explain the tail in the ρ_cal/ρ_obs distribution (Fig. 1) with low retrieved ρ_cal values.

Fig. 8. Scatterplots of hmF2 vs. LogNmF2 for all the analyzed IRO observations in a comparison with the climatological IRI-95 model for solar minimum and solar maximum conditions.

Figure 8 demonstrates one more inconsistency related to IRO observations. The IRO hmF2 vs. LogNmF2 dependence on average is a direct one, while IRI gives an inverse dependence: larger hmF2 correspond to lower logNmF2 values. The inverse type of hmF2 vs. NmF2 dependence in their diurnal variations is a correct one from physical point of view. It is due to the simultaneous influence of photo-ionization, neutral temperature, and thermospheric winds diurnal variations (e.g., Rishbeth & Garriott 1969).

All these results manifest the inconsistency between the COSMIC RO EDP and physics of the F2-layer formation which our method is based on. The absence of a solution tells us that there is no such a combination of aeronomic parameters which could provide the observed Ne(h) distribution. Therefore the proposed method can be used to check the physical correctness of the observed EDP. It should be stressed that observed IRO EDP may be absolutely correct from the point of view of the IRO method, but the majority of them are not usual Ne(h) vertical profiles used in ionospheric physics, and from this point of view such profiles should be considered as incorrect.

Continuing the analysis under solar maximum conditions, we will analyze how much the limitations in the altitude extent of CHAMP IRO EDP affect the performance of M2012 method.

In CHAMP IRO observations EDP are available below 400–410 km and the upper boundary condition is taken at 400 km height in our calculations. This is too close to the F2-layer peak and this limitation may affect the fitting process. It was supposed that poor testing results obtained for solar maximum may be related to the problem with the height of CHAMP satellite. To check this we took Millstone Hill ISR observations for October 12, 2002 and retrieved thermospheric parameters with the upper boundary conditions at 400 km and at 550 km. CHAMP neutral gas density observations are also available in the vicinity of Millstone Hill for this date at 18.08 UT.

Both runs gave good solutions with perfect Ne(h) fitting, but the retrieved parameters turned out to be different. The solution with the upper boundary at 400 km gave higher Tex = 1404 K and larger gas density ρ = 1.22 × 10⁻¹⁴ g cm⁻³ compared to Tex = 1362 and ρ = 1.02 × 10⁻¹⁴ g cm⁻³ obtained with the upper boundary specified at 550 km. It is interesting to note that this ρ = 1.02 × 10⁻¹⁴ g cm⁻³ exactly coincides with CHAMP neutral density observations reduced to the Millstone Hill location. Our analysis of all calculations for solar maximum gave average ρ_cal/ρ_obs = 1.20 (see earlier). This also coincides with ρ₄₀₀/ρ₅₅₀ ≈ 1.2 obtained in two calculations. So an average 20% bias obtained in our calculations for solar maximum can be related to the upper boundary condition specified at 400 km height. Therefore, CHAMP RO profiles available below 400 km cannot be used with our method to retrieve thermospheric parameters under solar maximum conditions. Additional IRO observations up to ≥600 km under solar maximum are required to check the efficiency of our method under high solar activity.

5. Discussion and conclusions

We are developing a new method to monitor the upper atmosphere neutral composition and temperature solving an inverse problem of aeronomy. It is well known that the observed EDP at F2-layer heights contains the whole information on thermospheric parameters, solar EUV, neutral winds, and electric fields and all these aeronomic parameters in principle can be extracted from EDP. The idea of this approach is very straight and clear but its implementation is not simple. An analysis of this problem may be found in Mikhailov & Lilensten (2004). Such a method has been developed and validated (Mikhailov et al. 2012, 2013) and its applicability in an operational environment has been assessed using various sources of EDP observations as input information. Tests were made using EDP from ISR, from Digisondes (Mikhailov et al. 2012; 2013), and in this publication we are using EDP from IRO. Today the main problem is seen in the accuracy of the input EDP rather than in the method itself. Wide testing of the method has shown that it demonstrates stable performance and gives acceptable results provided that accurate EDP are used as input. Incoherent Scatter Radar EDP measurements were shown to provide the necessary accuracy of observations both at middle latitudes (Millstone Hill ISR) and at the geomagnetic equator (Jicamarca ISR) while with autoscaled Digisonde EDP the method’ s performance is stable but less accurate comparing to ISR EDP.

Radio occultation is a relatively new and efficient method to investigate the Earth’s ionosphere and millions of EDP have been obtained for the past decade. These observations are widely used both for global ionosphere modeling and analysis of physical processes in the upper atmosphere. The possibility to ingest IRO EDP to the proposed method seems interesting as it may open wide opportunities for the thermosphere investigations.

The undertaken analysis has shown that middle latitude daytime IRO EDP technically can be used with the proposed method, under the condition that the hmF2 and NmF2 parameters extracted by the specific EDP are compatible with the mid-latitude F2 layer formalism adopted by the M2012 method. A comparison of the retrieved neutral gas density with CHAMP observations under solar minimum conditions, has shown that when the hmF2 extracted from the IRO profile is within the range of IRI-hmF2 values for such conditions, the inaccuracy that is obtained for the retrieved neutral density is close to the announced absolute uncertainty ±(10–15%) of the neutral gas density observations with the CHAMP satellite (Bruinsma et al. 2004). In parallel, for the same set of cases the JB-2008 model gives MRD = 16.4% and the MSISE-00 model gives MRD = 28.9%. This demonstrates that the method can provide successful results using IRO EDP.

However our analysis gave a large percentage (~50%) of EDP for which either none or incorrect solutions were obtained. Lee et al. (2012) have also reported large percentage (35%) of rejected FORMOSAT-3/COSMIC RO profiles.

According to Schreiner et al. (2007) the horizontally inhomogeneous irregularities in the neutral atmosphere and the small-scale irregularities in the ionosphere (whose effect is not completely eliminated by the ionospheric calibration) affect the accuracy and the precision estimated from pairs of closely collocated occultations when the separation of the occultation planes is larger than the correlation radii of the irregularities. In addition, horizontal variations in the ionosphere and local spacecraft multipath can also contribute to differences and may be the cause of the large-scale fluctuation seen near 400 km altitude.

These effects may explain the uncertainties identified in the relationship between NmF2 and hmF2, which cannot be predicted using the physical theory and which contradicts to the middle latitude F2 layer formalism, and therefore for these cases no solution can be obtained. A future step would be to provide some quantitative criteria that will allow operators to assess in real-time the quality of the IRO EDP. For instance the development of a statistical algorithm to check if the extracted hmF2 and NmF2 parameters are within the range of IRI predictions for the specific geophysical conditions, could provide such a metric, but this is a topic for a future investigation.

The main conclusions are summarized as follows:

1. The M2012 method that retrieves the main thermospheric parameters (Tn, [O], [N₂], [O₂]) from EDPs in the daytime mid-latitude F2-region was applied for the first time to IRO Ne(h) observations conducted under solar minimum (2007–2008) and maximum (2002) conditions.

2. It was shown that daytime mid-latitude IRO profiles technically can be used to retrieve the main thermospheric parameters. A comparison with CHAMP neutral gas density observations under solar minimum conditions has shown that, under conditions where the hmF2 parameter extracted from the IRO EDP is within the range of IRI predictions for the given geophysical conditions, the neutral gas density can be retrieved with an inaccuracy close to the absolute inaccuracy ±(10–15)% of CHAMP observations. For the same set of cases, empirical models JB-2008 and MSISE-00 being compared to the same observations provide less accurate results.

3. However IRO observations manifest a large percentage (~50%) of rejected EDP which cannot be used with our method due to their insufficient accuracy. The method is very sensitive to the quality of the EDP. Comparisons with IRI predictions and with coincident ISR EDP, indicate that inconsistency between NmF2 and hmF2 or incorrectness in the Ne(h) shape (inevitable in IRO observations) prevents the M2012 from providing a solution. The large percentage of rejections indicate that IRO profiles cannot be considered as any other type of Ne(h) vertical profiles and the mid-latitude daytime F2-layer formalism adopted by the M2012 method cannot be routinely applied to such profiles.

4. The performance of the method with IRO EDP under solar maximum conditions is in general more stable, and this matches well with the comparison of the hmF2 parameter extracted from the IRO EDP with the IRI predictions, where a much better agreement is seen. However, under such conditions the only available IRO EDP were observed from CHAMP at heights ≤400 km. This is not sufficient for getting a correct solution as the upper boundary should be specified at heights ≥550 km under solar maximum. This resulted in retrieved neutral gas densities with a 20% bias with respect to CHAMP observations, presumably due to incorrect specification of the upper boundary conditions.

5. A future step should be done to specify the quantitative criteria that will allow the automatic assessment of the applicability of the IRO EDP, most probably based on the IRI predictions. This will be a definite step for the implementation of the method in an operational environment, keeping in mind more IRO EDP available from future missions (e.g., COSMIC-2).

Acknowledgments

The authors thank the Millstone Hill Group of the Massachusetts Institute of Technology, Westford for running the radar and providing the data. We also acknowledge NSPO of Taiwan and UCAR of USA for providing the FORMOSAT-3/COSMIC data and the ISDC of GFZ for the CHAMP data. The research was supported by the NATO Collaborative Linkage Grant CLG 984141 and the COST Action ES0803.

References

Appendix

All Tables

All Figures

	Fig. 1. Distributions of R = ρ cal/ρ obs ratio for the retrieved cases with the M2012 method (top panel), and the corresponding distribution based on JB-2008 and MSISE00 neutral gas densities models for solar minimum conditions (middle and bottom panel). The average R ave and the standard deviation SD values are given. The area between dashed lines corresponds to R ave ± SD.
In the text

	Fig. 2. Retrieved and modeled neutral gas densities from Table A1 are given vs. the observed values. Note that the retrieved densities are more centered with respect to the observed ones while modeled values are biased overestimating the observations. MRD and bias are given for a comparison.
In the text

	Fig. 3. Same as Figure 1 but for solar maximum conditions.
In the text

	Fig. 4. Scatterplots of hmF2 and NmF2 determined by Millstone Hill ISR vs. RO method.
In the text

	Fig. 5. Daytime Millstone Hill ISR EDP (triangles) in a comparison to COSMIC RO measurements (solid line) under solar minimum 2007–2008 conditions. Six of the best (top panels) and six of the worst (lower panels) coincidences are shown.
In the text

	Fig. 6. Observed with COSMIC and Millstone Hill ISR EDP on January 31, 2008.
In the text

	Fig. 7. Observed at Millstone Hill ISR in July 2008 daytime hmF2 vs. LogNmF2 in a comparison with COSMIC cases for which the solutions are within the R ave ± SD range and cases for which no solutions were obtained.
In the text

	Fig. 8. Scatterplots of hmF2 vs. LogNmF2 for all the analyzed IRO observations in a comparison with the climatological IRI-95 model for solar minimum and solar maximum conditions.
In the text