© A. Belehaki et al., Published by EDP Sciences 2012

1. Introduction

The safety and security of space operations requires that spacecraft operators know and understand the environment around their spacecraft. A measure of the ionospheric element of that environment is given by the precise values of the total electron content (TEC) between the spacecraft and some location on the surface of the Earth (partial TEC) – and in some cases between two spacecrafts. The partial TEC value differs from those derived from ground-based GNSS receivers. The latter values represent the TEC of the whole ionosphere-plasmasphere system, i.e. height-integrated from the surface to the GNSS satellites around 20 000 km altitude. For many space activities we must partition GNSS TEC above and below the spacecraft in order to derive the TEC applicable to a particular activity. A precise knowledge of the electron density distribution with height is required to calculate the TEC between two arbitrary points in space.

As direct observations of electron density are sparse, precise estimation of the partial TEC requires accurate modeling of the topside ionosphere, i.e. the region above peak electron density and that gradually merges into the plasmasphere at altitudes above the O⁺-H⁺ transition height. The topside is the region that contains most of the plasma that contributes to the TEC values that we seek to estimate; typically only 25% is below the F2 peak while typically <10% comes from the plasmasphere. Furthermore, the topside ionosphere is the region in which most LEO spacecraft operate, so accurate estimates of the TEC relevant to a particular spacecraft activity require us to partition the TEC within this region. Thus modeling the topside electron density distribution is central to our key objective.

Along this direction, Marinov et al. (2004) and subsequently Kutiev et al. (2006), Kutiev & Marinov (2007), and Belehaki et al. (2006a) developed an improved topside profiler model, the Topside Sounders Model assisted by Digisonde (TaD). The first version of TaD (TaDv1) is based on the Topside Sounder Model (TSM) which is an empirical model for the O⁺-H⁺ transition height (h_T), the topside electron density scale height (H_T), and their ratio R_T = H_T/h_T. TSM equations are based on vertical electron density profiles (EDP) measured by the topside sounders onboard Alouette-1a, 1b, 1c, and 2, and ISIS-1 and -2 (Bilitza 2001). The two topside parameters were extracted from 176,622 measured EDPs by assuming that the lowest gradient in EDP represented the O⁺ scale height, while the transition height was obtained by extrapolating the O⁺ density profile upward to the height where its density became one half of the electron density. The models provided the topside scale height H_T and O⁺-H⁺ transition height h_T on a global scale, as a function of geomagnetic latitude, day of year, local time, solar flux F10.7, and Kp. Kutiev & Marinov (2007) have found that O⁺ scale height and transition height in individual profiles correlate well, as in midlatitudes correlation exceeds 0.8. That property has been further used to develop a new model TSM, combining the former H_T and h_T models, with the model of their ratio R_T = H_T/h_T. The model ratio R_T, which is based on the individual ratios of scale height and transition height, can be used as a link to other applications. Kutiev & Marinov (2007) have also proposed a new profiler (TSMP), which used TSM topside parameters to provide the shape of the Ne profiles above the F2 layer peak height (hmF2). The TSMP (TSM profiler) describes EDP as a sum of O⁺ and H⁺ profiles. In the topside F region, the O⁺ profile is presented by the α-Chapman formula with a scale height H_T, and above the transition height the H⁺ distribution is presented by a pure exponential function of height with a scale height Hp = 16H_T.

Kutiev et al. (2009a) applied TSMP to Digisonde measurements to reconstruct the electron density profile from the F layer peak up to GNSS orbits developing the TaD model based on a new expression for the H⁺ scale height in the plasmasphere (Hp) extracted from ISIS-1 topside sounder data, as a function of the geomagnetic latitude. TaD method is considered to be a powerful basis for the development of operational applications especially over regions covered by dense Digisonde networks, as, for example, the European middle latitude region covered by DIAS network (Belehaki et al. 2006b; Belehaki et al. 2007), but in principle, the successful transition from research to operational models requires among others the systematic assessment of the models’ performance under all possible conditions (Tsagouri et al. 2012). The evaluation of the TaD model’s performance has been previously attempted in qualitative terms by Belehaki et al. (2009). The tests were mainly focused on the comparisons between the integral of the TaD profiles and the TEC measured by the CHAMP satellite (from 400 km up to GNSS orbits) and by ground-based GPS receivers. Comparison results of few IMAGE/RPI plasmagrams with TaD modeled profiles from conjugated Digisondes with the satellite have also been considered. Although a general agreement between TaD results and the corresponding parameters derived from the processing of satellite observations (CHAMP and RPI/IMAGE) and ground-based GPS receivers data was found, one has to note that this is just a qualitative indication of the model performance, since all these are comparisons made with indirect measurements of plasmaspheric parameters (electron density and TEC) and were mostly restricted to case studies. Therefore, a comprehensive validation plan based on the systematic comparison between TaD-derived and measured electron density profiles, especially in the topside ionosphere, is still a strong requirement.

In addition, Kutiev et al. (2006) have shown that the range of altitudes where O⁺ distribution is pure exponential is 165 km, while Kutiev et al. (2009b) showed that the H⁺ distribution is practically exponential above 2500 km, the accuracy of the TSMP profiling approach in the altitude region between the two exponential parts of the Ne profile is not yet extensively investigated. In particular, it is the assumption that H⁺ is the major light ion throughout the region of interest that still needs to be verified upon the measured profiles.

Starting from the latter, the investigation of the validity of the underlying TSMP assumptions in the transition region between topside F region and plasmasphere is the first point to be addressed in this paper. The results of this investigation guided a substantial reformulation of the underlying TSMP expressions, and the optimization of the TaD algorithm. Improvements of the TaD prediction capabilities in respect to the previous version of the model are verified through a self-test procedure based on the comparison between TaD and ISIS-1 topside profiles. Its accuracy at the topside is investigated through the comparison between TaD predictions with actual profiles obtained from Malvern-ISR.

2. Upgrade of the TaD model

2.1. Electron density distribution in transition region between topside F region and plasmasphere

The following discussion is based on the analysis of 12 616 Ne profiles measured by the topside sounder onboard ISIS-1 satellite during the period 1969–1971. The database was compiled after selecting those profiles, for which the scale heights H_T and Hp, and the transition height h_T lie within predefine limits.

As an indicative example of measured and calculated profiles, Figure 1 shows one measured Ne profile (red curve), the inferred O⁺ distribution (light gray line), and the transition height (yellow line). In principle, TSMP determines the H⁺ distribution (dark gray line) as having its maximum value at h_T equal to that of O⁺ and scale height Hp is obtained from the upper part of Ne. It is clearly seen in this example that the H⁺ values systematically exceed the measured Ne ones at altitudes above 1500 km. If one draws the H⁺ distribution at its right position, i.e. over the Ne profile (dashed green line), its maximum value at the transition height is significantly less than that of O⁺. Obviously, there is an excess of ionization, which complements the H⁺ density to become one half of the measured Ne. We speculate that this excessive ionization is assigned to He⁺, which is the only light ion that is expected to be found at these altitudes, besides the H⁺. Indeed, the behavior of the light ion density in the topside F region and the plasmasphere has been extensively studied by satellite data and in general, the main physics is known. Since Moffett & Hanson (1973) have proposed excessive He⁺ ions at low latitudes as formed through the E × B drift at equator, several basic studies have been conducted to reveal He⁺ distribution, among them Heelis et al. (1990), Kutiev & Stankov (1994), and recently Su et al. (2005). The blue curve represents the expected light-ion (H⁺ + He⁺) distribution.

Fig. 1. One measured Ne profile (red curve), the inferred O+ distribution (light gray line), the transition height (yellow line), and the H+ distribution above transition height (green dashed line). The TSMP H+ distribution is plotted with the dark gray line. The blue curve represents the expected light-ion (H+ + He+) distribution.

To further investigate the above finding, Figure 2a presents three daytime Ne profiles measured at three different geomagnetic latitudes: one in the equatorial zone (2° N), one in the lower middle latitude zone (43° N), and one in the higher middle latitude zone (50° N). The other three panels of Figure 2b, c, and d) analyze each of the profiles and the inferred ion distributions in the same format as in Figure 1. Based on Figure 2, the profiles in the equatorial and the higher middle latitude zones seem purely composed by H⁺ and O⁺, providing no evidence for the existence of extra ionization at these latitudes. However, the extra-ionization signature is clearly seen in the lower middle latitudes, indicating that the appearance of He⁺ is latitude dependent and most probably observed at middle latitudes.

Fig. 2. (a) Daytime Ne profiles measured at three different geomagnetic latitudes: equatorial (green), midlatitude (red), and higher midlatitude (blue). (b, c, and d) The figures in each of the three panels reproduce each of the profiles and the inferred ion distributions in the same format as in Figure 1. The two purple arrows in panel c formulate the difference between the inferred H+ and O+ densities at the transition height as evidence for the existence of extra ionization. Single purple arrows in panels b and d indicate coincidence between the two values.

Next, all ISIS-1 profiles that were included in the Alouette/ISIS database were analyzed in a similar manner as above approach. In particular, the difference ln(O⁺) − ln(H⁺) at the transition height was estimated for all cases and the results are presented in Figure 3a as a function of the geomagnetic latitude. The difference, which here is attributed to the presence of He⁺, is larger at middle to low latitudes and less pronounced at higher latitudes. Another view of the He⁺ presence is given in Figure 3b, which provides the geomagnetic latitude/local time contour map of the difference between O⁺ and H⁺ densities at the transition height. The difference is assumed to be directly equal to He⁺ density. This map clearly shows that at middle latitudes the He⁺ appears in the morning and the afternoon hours, while it is almost absent at noon and at night. Contrary, at high latitudes He⁺ appears mainly at night.

Fig. 3. The difference ln(O+) − ln(H+), for all ISIS-1 profiles, taken at the transition height as a function of geomagnetic latitude is shown on the left side. The geomagnetic latitude/local time contour map of the difference of O+ and H+ densities at the transition height is given at the right panel.

There could be other explanations of the inconsistent exponential distribution in transition region, for example, a vertical electrodynamic drift at equator, which increases O⁺ in the anomaly crests and hence makes its vertical distribution to deviate from exponential. There is no direct way to verify which assumption is correct: the existence of He ions or the effect of equatorial E × B drift which alters the O⁺ distribution around the ionospheric anomaly crests. However, if the vertical electrodynamic drift at the equator changes the O⁺ distribution at the crests, it will necessarily change the He⁺ distribution creating a latitude maximum poleward from that of the O⁺ crests. Therefore, the drift is expected to increase He⁺ density in the same proportion as it does for O⁺, but at little higher latitudes. The net effect is again the increase of He⁺ density.

2.2. Improvements in the TSMP formulation

Based on the above discussion, one may argue that it is possible to accommodate the observed extra ionization due to the presence of He⁺ ions and its characteristics in the TSMP formulation. For this purpose, a correction parameter g, defined as the ratio of the H⁺ and O⁺ densities, g = n(H⁺)/n(O⁺) and a new term similar to the second term describing the H⁺ distribution were introduced in the original TSMP formula (e.g. Kutiev et al. 2009b) that provides the electron density as a function of the altitude. The new formula is now given in the form:(1)

In the absence of He⁺, g = 1 and the third term is equal to zero. The lower the g is, the more He⁺ ions are present. Obviously, g should be a function of the geomagnetic latitude and the local time and probably the season.

Concerning g, TaD cannot obtain any information about the He⁺ from the bottomside sounding results and therefore its presence is a priori determined on the basis of the topside sounder profiles. To obtain analytical expression for g, profiles from the ISIS-1 database were used as it is described previously but excluding those with g = 1. This restriction is necessary in order to exclude contribution of the profiles with no presence of He⁺. To approximate g obtained from the individual measured profiles (g = n(H⁺)/n(O⁺)), we chose analytical function of three variables: local time, geomagnetic latitude, and O⁺ density. The latter variable was added to the first two after finding the strong correlation between g and n(O⁺). A description of the mathematical formulation for g approximation is provided in the Appendix A.

The scatter plot of model g (gm) versus measured g values is given in Figure 4. The thick blue line is the linear regression through the origin over gm values. Correlation coefficient is 0.98, which shows that the analytical function of g represents the data exceptionally well. However, the angular coefficient of the regression line is 0.945, which means that at higher g values gm tends to overestimate the g function. The parameter g is a 3-D function that could not be easily visualized. To get an idea of how g function depends on each of the three variables, in Figure 5 we have approximated the projection of gm values on each of the axes plane with low order polynomials. In the top panel, the linear fit shows that gm tends to decrease with increasing O⁺ density. This means that He⁺ (reversely proportional to g) predominantly appears when O⁺ density increases, which happens at the F region anomaly crests. The middle panel shows that gm has a minimum of around 0.6 at low latitudes and increases toward higher latitudes, which complies with the known fact that He⁺ is predominantly present at low latitudes (Fig. 3, left panel). Local time dependence shows a minimum in the early afternoon hours, when the equatorial anomaly is well developed. As a conclusion, we may say that gm behavior is in compliance with expected variation of He⁺ density and the function can reliably estimate contribution of helium ions in the transition region.

Fig. 4. The scatter plot of the modeled g (gm) versus measured g values. The thick blue line is the linear regression through the origin over gm values.

Fig. 5. The projection of the function gm on each of the axes plane approximated with low order polynomials; top panel: the linear fit shows that gm have a tendency to decrease with increasing O+ density; middle panel: gm has a minimum of around 0.6 at low latitudes and increases toward higher latitudes; bottom panel: local time dependence shows a minimum in the early afternoon hours.

2.3. Optimization of the fitting procedure

Regarding the scale height of He⁺, it is reasonable to assume that He⁺ has its maximum again at the transition height, h_T, and exhibits a pure exponential distribution above, with a scale height equal to four times that of H_T. However at higher latitudes, for the plasmaspheric scale height Hp, which theoretically should be 16 times that of H_T, the multiplication factor drops significantly due to the shape of plasmasphere (Kutiev et al. 2009b). Having in mind this fact, we tried to infer the factor from the data. So, formula (1) is modified as in the last exponent instead of 4H _T we place sH _T:(2)

By definition, the He⁺ density strongly depends on the O⁺ and H⁺ densities at the transition height, which in their turn strongly depend on the scale heights H_T and Hp, determined from the lowest and the highest altitude gradients of the measured Ne profiles. Increasing the fitting accuracy of the measured Ne profiles in the transition region comes up through the increase of the accuracy in the determination of the He⁺ at h_T. To increase the accuracy of the Hp and H_T, we developed a sophisticated algorithm, based on the optimization of H_T, Hp, and s toward best fitting with measured Ne profile. Formula (2) runs multiple times by using a set of H_T, Hp, and s and the model profile is compared with the measured Ne profile by calculating the root mean square (RMS) error in each run. A series of values is assigned for each of these variables and Ne is calculated from formula (2), executing three nested cycles: the outer for H_T, middle for Hp, and internal for s. Each run takes one of the H_T values, calculates the corresponding H_T, takes one of the Hp values, calculates the ratio g at h_T, takes one of the s values, and then calculates the whole profile based on formula (2) in the altitude range corresponding to the measured profile. The RMS error is stored and the next cycle is performed. We have set 40 values around the minimum H_T, 40 values around maximum Hp, and 120 values of s around 4. The whole optimization of a given profile contains 40 × 40 × 120 = 192 000 runs. So, we select the set of the three variables that has the smallest RMS error. In this way we perform optimization to 14 775 measured ISIS-1 profiles.

We found that the He⁺ scale height cannot be reliably obtained. Indeed, the RMS error is formed by comparing formula (2) with the whole measured profile. He⁺ density in some cases is comparable with electron density within the transition region, but it is always much smaller outside that region and its net contribution to the RMS error is negligible. When calculating procedure executes formula (2) for variables in the nested cycles, RMS errors are practically determined by the pairs of H_T and Hp values, with no significance of values of s. Thus s value could not be optimized. Statistically for the whole database, s values are spread almost uniformly in the defined range. Taking into account this fact, we further take s = 4.

It is well accepted (see, for example, Banks & Kockards 1973; Stankov et al. 2003) that the broad maximum of the vertical H⁺ distribution is placed slightly above the transition O⁺-H⁺ height, e.g. above the transition region with the highest rate of charge-exchange reaction O↔H⁺. The maximum of the vertical distribution of He⁺ is found to lie around or slightly below transition height (Heelis et al. 1990; Koleva & Kutiev 1990; Kutiev & Stankov 1994). To simplify the TSMP model, we placed both H⁺ and He⁺ maxima at the transition height. We assumed that should not worsen the accuracy of the profiler because TSMP reconstructs the total Ne profile and the ion profiles only compose it. TSMP is not designed to provide detailed ion composition from the Ne profile.

Figure 6 shows, as an example, the profile measured on day 97 of the year 1969, at 17:03 UT around equator (lat = 12° S, long = 12.5° E). The ion density profiles in logarithmic scale obtained without optimization (denoted as “old”) and those with optimization (denoted by “new”) are color coded, along with the measured Ne profile presented with the black line. The “old” profiles are presented by dashed lines, while the “new” profiles are given by solid lines. The total ion density presented the sum of the “new” profiles is given by the orange line. In this sample, optimization increases slightly the values of H_T and h_T and reduces Hp. Optimization invokes its largest effect on He⁺ density, which is substantially reduced. However, He⁺ density profile below transition height h_T depends strongly on H⁺ density shape, which, as we have demonstrated, could not represent the real situation. As we are interested in reconstructing the electron density profile, which below h_T is mainly defined by O⁺, the shape of H⁺ is not important. The same is true for the He⁺ profile below h_T. In this study we are mainly interested in He⁺ density around its maximum fixed at h_T. The standard deviation of the optimized model from measured profiles, denoted as RMS, is 0.051, while without optimization it is 0.112 units of logarithm density. Here optimization has reduced RMS error more than twice.

Fig. 6. A sample profile measured on day 97 of the year 1969, at 17:03 UT around equator (lat = 12° S, long = 12.5° E). The ion density profiles in logarithmic scale obtained without optimization (denoted as “old”) and those with optimization (denoted by “new”) are color coded, along with the measured Ne profile presented with the black line. The “old” profiles are presented by dashed lines, while the “new” profiles are given by solid lines. The total ion density presented the sum of the “new” profiles is given by the orange line.

3. Verification of the model’s performance

Verification aims at determining how well an operational model is performing. To verify the performance of the improved TaD, we compare results of the improved TaD (TaDv2) with actual measurements obtained by ISIS-1 satellite, and we compare the RMS errors obtained with the previous (TaDv1) and the new version (TaDv2) of the model to confirm improvement in the performance.

This test, as it is based on comparison with ISIS-1 measurements, is not an independent evaluation test, because ISIS-1 profiles were used for the development of the TaDv2, but provides with evidence of how well the new profiler reproduces the input measurements. It is also designed to verify improvements in respect to the previous version of the TaD model (TaDv1).

The comparisons were made through the estimation of the RMSE over whole profile, according to the following formula:(3)where np indicates the number of points in each ISIS1 EDP, h_min and h_max are the lowest and highest altitudes for which ISIS1 electron density measurements are provided for a given profile, n_TaD and n_ISIS1 are the electron density values at a given height h, provided by the TAD model and ISIS1 satellite respectively.

The distribution of the RMS errors comparing modeled with measured profiles before (no optimization) and after (with optimization) TaD upgrade collected from the whole selected ISIS-1 database is presented in Figure 7 in boxplot and histogram formats. The boxplot format (left panel) includes a box and whisker plot for each case. The box has lines at the lower quartile, median (red line), and upper quartile values. Whiskers extend from each end of the box to the adjacent values in the data; in our case to the most extreme values within 1.5 times the interquartile range from the ends of the box.

Fig. 7. The distribution of the RMS errors comparing modeled with measured profiles before and after TaD upgrade for the whole selected ISIS-1 profiles in: (a) boxplot format (left panel) and (b) histogram format (right panel).

The histogram distribution of the RMS errors is presented in the right panel. The green filled bars represent the histogram of the RMS error for the previous TaD version, while transparent blue bars show the distribution of RMS error for the improved profiler. It is clear that the introduced upgrades improve both the accuracy and the precision of the model’s predictions. The average RMS error from 0.189 is now 0.085, being reduced 2.25 times.

An important evaluation test applies to the TEC parameter. For this purpose, TaD and ISIS1-derived TEC values were compared in two altitude zones: in the topside ionosphere and in the plasmasphere region. TEC was calculated by integrating the electron density from the lower height of the ISIS1 satellite (slightly above hmF2) up to the transition height and from the transition height up to the maximum height of the ISIS1 satellite, respectively, for the two zones. This test was applied to both the original and the improved TaD profilers.

The results are presented in Figure 8 in a scatter plot format. As a first comment, we can stress out the very good agreement between the TEC parameters estimated by the ISIS1 profiles and those provided by the improved TaD model reflected in a correlation coefficient greater than 99% for both altitudinal zones. Figure 8 also demonstrates obvious improvement in the model’s performance in respect to its previous version.

Fig. 8. Scatter plots of the TEC estimates derived from both ISIS-1 and TaD profiles, separately for the previous (left) and the improved (right) TaD versions and two altitude zones, the plasmaspheric (top) and the topside ionosphere (bottom).

To further investigate the relative performance of the two TaD versions, the distributions of the error between the model and the data-derived TEC estimates for all of the cases are presented in Figure 9 in a boxplot format. The results indicate substantial improvement in both the accuracy and the precision of the model’s predictions in TaD-TEC estimates for the plasmaspheric region. In the topside ionosphere, the improved TaD does not yield in noticeable improvement to the main features of the error distribution, but it succeeds to significantly reduce the number and the range of the outliers denoted with the gray crosses in the figure. This is expected since TaD improvements reported here apply mainly in the region above the transition height, so we expect that the performance of the previous and current versions of the model at the topside is equivalent. The He⁺ profile introduced in the TaDv3 becomes a significant component of the total density profile in the heights above the transition height and this is reflected in the improved performance of TaDv3 demonstrated with the results of Figure 9.

Fig. 9. Error distributions between the model and the data-derived TEC estimates for all of the cases considered in the paper.

4. Validation of the model in the topside part of the profile

Validation aims at determining that the model is ready for operational use, based on its performance. The topside ionosphere is a region of special interest because it is the region where inhomogeneities and scintillations occur predominantly, with the primary disturbance region being in the F region between 250 and 400 km. In addition, a large number of LEO satellites operate in that region and reliable information on the existence of density gradients is very important especially during disturbed space weather conditions. Therefore it is important to assess the performance of the model in this specific region. This is achieved with the comparison of the reconstructed EDPs with independent measurements of the electron density distribution from the ISR operated in Malvern in 1968–1971. This data set has been chosen because Malvern is a middle latitude site and the operation of the ISR temporary coincides with the ISIS-1 mission, therefore long-term trends in the ionosphere cannot be the cause of any possible discrepancies in the comparison results. In addition, validation of space weather operational models with ISR profiles has been recommended by the US National Space Weather Program for the establishment of the top high priority metric for the ionosphere-thermosphere domain, as being the only independent source of full density observations from the ground.

It has to be emphasized here that this test does not mean to assess improvements made in TaDv2 regarding the contribution of the He ions tn the reconstructed profile, since this concerns the part of the profile above the topside. Based on the conclusion of the previous section in the topside, both versions of the model have an almost equivalent performance, and therefore we present here results only for TaDv2. This test aims at providing information on the TaD performance in quantitative terms, which has not been done before and this information is very important for users willing to apply the model in the altitude zone of 200–600 km (approximately).

For this purpose, a data set of some of 4000 electron density profiles obtained at Malvern-ISR site (52.1° N, 2.3° E) during the time interval 1968–1971 were analysed. A quality check was applied to the data, to separate profiles with missing topside part and those with large scatter in its topside part which does not allow determining reliably the scale height. After the quality check, almost the 30% of profiles have been eliminated from our sample. The F layer peak density NmF2 and height hmF2 were extracted from each individual profile. However, for the scaling of the ISR EDPs we have applied the method used in the topside sounder profiles scaling, i.e. the NmF2 and hmF2 parameters have been located in each measured profile; then its topside scale height H_TM has been obtained as a regression line from a number of data points above the peak, with the transition height h_TM being calculated from H_TM and corresponding TSM ratio R_T. The four quantities (foF2, hmF2, H_TM, and h_TM) taken from the measured profile have been ingested to the formula (2) and the model profile was calculated for altitude range of the topside part of measured profile. Keeping this in mind, the analysis attempted below may be considered as a rough evaluation of the method’s performance but in our opinion it provides suggestive conclusions. Figure 10 gives some examples of the derived profiles.

Fig. 10. Examples for TaD-derived profiles based at Malvern site. The ISR EDP are denoted with the cross symbol. Red line indicates the modeled O+ profile and blue line indicates the modeled EDP.

Malvern-ISR profiles cover the altitude range up to approx. 700 km which is about the altitude of the actual transition height, so it is expected that O⁺ will dominate and the observed Ne profile will follow closely the TaD profile of the O⁺ density. Based on this, we first compared the model predictions for O⁺ distribution in respect to the measured O⁺ profiles. The results are expressed in terms of the normalized RMSE estimates that were obtained over the topside profiles (from hmF2 up to the maximum height of the observed profiles). Their distribution is presented in Figure 11.

Fig. 11. The normalized RMSE distribution of model predictions for O+ in respect to the measured profiles.

The distribution of normalized RMSE for the predictions of the total electron density in respect to the measured TEC from Malvern-ISR is given in Figure 12. Comparing the results in Figures 11 and 12 we note a slightly better fit between the observed and modeled O⁺ densities rather than observed and modeled total electron densities. This is probably due to the fact that the upper part of ISR profiles is not included in calculation of H_TM, because as explained above, the topside scale height H_TM has been obtained as a regression line from only a number of data points above the peak. In addition, due to the fact that the topside ISR profiles are limited to 700 km, transition height cannot be always reliably obtained, e.g. we cannot reach the height where O⁺ is one half of the measured density. Moreover, the measured profiles around their upper height values exhibit often an artificial increase which confuses the TaD algorithm.

Fig. 12. The normalized RMSE distribution of modeled total electron density in respect to the measured profiles.

This might be the reason for the tendency of the model to overestimate the electron density, as seen in Figure 13. Here we present the distribution of the simple difference between the observed and the modeled electron densities at three selected altitudes: 400 km (top), 500 km (middle), and 600 km (bottom). Another interesting result concerning the dependence of the model’s performance with the height is that the fit of the model to the observed values increases with increased altitude. This can be explained taking into account that the density decreases with altitude and thus the simple difference becomes smaller.

Fig. 13. The distribution of the simple difference between the observed ISR and the modeled electron densities at 400 km (top), 500 km (middle), and 600 km (bottom).

As TEC is the most critical parameter for the reliable performance of several applications, we also performed comparison tests between TEC estimates obtained from modeled and observed EDP from the ISR at Malvern site concentrating on the topside part of the profile, i.e. above the altitude of hmF2.

First, we present in Figure 14 the scatter plot of TEC modeled versus TEC observed estimates. The results indicate that the two parameters correlate reasonably well. A larger scatter is seen in the daytime hours (the area with larger TEC values), which is probably due to the high altitude of the transition height, something that makes its accurate determination problematic due to the lack of observations above 700 km. The O⁺/H⁺ transition height varies but seldom drops below 500 km at night or 800 km in the daytime, although it may lie as high as 1100 km, depending on the geophysical conditions, and particularly on solar activity (Denton et al. 1999).

Fig. 14. The scatter plot of TEC modeled versus TEC observed estimates.

In Figure 15 we present the TEC error (abs(TECobs-TECmod)) distribution and the relative TEC error distribution. At the top we present the results based on all available measurements, in the middle we present the corresponding distribution diagrams for a subset of measurements that correspond to daytime EDP, and in the bottom we give the distribution diagrams for nighttime cases. The analysis on the whole sample of measurements gives a mean error of 3TECU and mean relative error of 28.5%. However, this error is considerably reduced when we include only the nighttime cases. This is an additional evidence showing that during daytime, when the transition height is in its upper limits, or even higher of 700 km, there are not enough data for the TaD to estimate correctly the transition height and this yields to artificial results. The same problem was reported by Belehaki & Kersley (2006), when they compared the ISR EDP from Malvern with the results from the topside extrapolation based on the Reinisch & Huang model (2001).

Fig. 15. The TEC error (abs(TECobs-TECmod)) distribution and the relative TEC error distribution. From the top to the bottom: (a) results based on all available measurements; (b) the corresponding distribution diagrams for a subset of measurements that correspond to daytime EDP; (c) the distribution diagrams for nighttime cases.

5. Discussion and conclusions

The performance of the first version of the TaD paper was evaluated in qualitative terms by Belehaki et al. (2009) through the comparisons between electron density profiles and TEC parameters extracted from TaD model and (a) CHAMP reconstructed profiles, (b) CHAMP-derived TEC parameters, (c) ground-based GPS-derived TEC parameters, and (d) profiles reconstructed from RPI/IMAGE plasmagrams. All data that have been used come from indirect measurement of the related parameters (electron density and TEC) and therefore, due to the uncertainties involved in all data sets it was possible to extract only results concerning the qualitative characteristics of TaD. We have reported a general agreement between TaD results and the corresponding parameters derived from the processing of satellite observations (CHAMP and RPI/IMAGE) and ground-based GPS receivers data. The preliminary analysis of the TaD performance during storm intervals indicated that TaD is extremely sensitive to the Digisonde measurements and to disturbances that reflect changes in the ionization in the bottomside ionosphere, indicating that TaD has the potential to provide qualitative characteristics of ionospheric-plasmaspheric conditions over Digisonde locations.

In this paper we propose further improvements to the model, attempting also a systematic validation that allows us to have a quantitative estimate of the model error. The improvements come from (i) the introduction of the He⁺ profile in the reconstruction of the total electron density profile, whose presence is more pronounced at the middle latitudes in the morning and afternoon hours and (ii) an approximation of the plasmaspheric scale height as a function of altitude, latitude, local time, and season, using an optimization procedure to best fit model with observations.

In general, He⁺ is a minor ion constituent, but in some cases its density is comparable with that of O⁺ in the transition region. He⁺ term in formula (2) allows a better fitting of the formula to the measured profiles, since it results in lower RMS. This is done through the optimization procedure. Inclusion of He⁺ term reduces the model error (increases accuracy) which is the main goal of the modeling. It also demonstrates the capability of the model to infer more information from the data and enlarge its range of applications. For example, if measured vertical Ne profiles in ionosphere and plasmasphere may become available in the future, TaDv2 model is capable to separate ion distributions and allow studies in geophysical sense.

The evaluation of the improved TaD performance verified significant improvement in respect to its original version in both the accuracy and the precision of the electron density and TEC estimates. The improved profiler reproduces with very high accuracy the ISIS-1 profiles demonstrating that the model does what it is designed to do and that it is free of intrinsic uncertainties. The very high (more that 99%) correlation coefficients between modeled and measured plasmaspheric and topside TEC is an additional confirmation of the model’s successful performance in the whole altitudinal range, from the topside ionosphere to the GNSS height.

We have also performed systematic comparison with 3 years of electron density profiles data, obtained from the Malvern-ISR. This comparison gave us the possibility to estimate the model error in the area of the topside ionosphere which is of particular interest for the operators of LEO and MEO satellites. The mean RMSE obtained from the comparison of ISR EDPs is in general larger than the RMSE obtained from the comparison with ISIS-1 profiles but this was expected as ISIS-1 profiles were also used for the model’s development. The comparison with ISR EDPs gives a model’s error mainly attributed to limitations imposed by the measured profiles themselves. During nighttime hours, when the transition height is lower than 700 km and can be determined by the TaD algorithm, the model values are in very high correlation with the ISR observations, with a mean error of only 1.1 TECU. An overall model error of 3TECU is estimated and it is close to the measurement (GNSS) error.

Further improvements on TaD model may be foreseen in two directions. First, based on recent developments giving evidence that the use of the variable scale height can simulate successfully electron density structures variation above the peak height (Reinisch et al. 2007; Kutiev et al. 2009b), one may argue that substantial improvement to TaD profiler may be achieved through its realization by using a variable scale height function instead of the TSMP approximation that uses the classical α-Chapman shape to represent the O⁺ distribution above the F2 peak hmF2 with a constant scale height provided by Digisonde measurements. For the description of F region profiles Rishbeth & Garriott (1969) suggested the use of a modified α-Chapman function, which takes account of the height variation of the scale height. Reinisch & Huang (2001) had used vary-Chap functions to represent measured bottomside profiles and determined from it the Chapman scale height Hm at the F2 peak. Height integration of the bottomside profile and the topside Chapman function provides an estimate of the ionospheric total electron content ITEC that is routinely calculated in all Digisondes (Reinisch et al. 2005). Of course by its derivation, ITEC does not include the plasmasphere electron content (Belehaki et al. 2003). Using the available plasmasphere and topside profile measurements from ISIS-1 satellite we tried to investigate the possibility to construct a suitable scale height function that changes with height for a topside vary-Chap function. Our analysis is described below and was based on the ISIS-1 database used also in the previous sections.

Figure 16a shows five Ne profiles in semi-logarithmic scale (ln(Ne), h) measured at lower middle latitudes in the afternoon hours. Arrows indicate the corresponding transition heights. In Figure 16b we present the derivatives dh/dln(Ne) which present the measured scale height versus altitude up to 1200 km and demonstrates how the scale height of the topside F region is obtained. Although the profiles themselves look pretty smoothed (Fig. 16a), obtaining the scale height is not so trivial technique, since the proper number of data points for linear regression should be selected automatically. In certain cases the technique fails to extract properly the scale heights and yield unreliable model parameters. That is why, it needs to be further tested and improved.

Fig. 16. (a) Five Ne profiles in semi-logarithmic scale (ln(Ne), h) measured at low-midlatitudes in afternoon hours. Arrows indicate the corresponding transition heights. (b) The augmented lower part, demonstrating the large variability of the scale height at the topside F region.

Figure 17 summarizes measured derivatives (measured values marked by red crosses) of the profiles from the whole database. The histogram of Hp obtained from the corresponding profiles is shown at the top portion of the figure. The average Hp value is around 1100 km and the scatter (standard deviation) is 40%. As was mentioned by Kutiev et al. (2009b), Hp is latitude dependent and this dependence is contained in the statistics on Figure 17. The large scatter of dh/ln(Ne) clearly indicates that an analytical expression of the varying scale height should contain dependences (besides altitude) on latitude, and probably local time and season. Without such a diversification, the total error (standard deviation of data around the model curve) will be too large and will compromise the reliability of reconstruction technique. We conclude that, at least on this stage, approximation of the altitude varying scale height with single expression is not feasible for TaD improvement due to the limited availability of data.

Fig. 17. A statistical representation of the measured derivatives (measured values marked by red crosses) of the profiles from the whole database. The histogram of Hp obtained from the corresponding profiles is shown at the top portion of the figure.

In another direction, studies performed during the last years demonstrated that the exploitation of coordinated measurements from ionospheric sounders and GNSS receivers can lead to the development of advanced algorithms for the accurate specification of the electron density distribution from the bottomside ionosphere to the plasmasphere. Following this current trend, further improvements for TaD model are attempted in an accompanied paper submitted in the current issue (Kutiev et al. 2012), where we propose adjustments of the integral of the modeled electron density to the TEC derived from GNSS signals.

Acknowledgments

This study has been supported by the COST Action ES0803. The authors also acknowledge the contribution of the EOARD through the Grant Award FA8655-12-1-0004. Anna Belehaki is grateful to Prof Len Kersley for providing the digitized electron density profiles from the Malvern-ISR.

Appendix A

Mathematical formulation

Coefficient matrix CFN is a solution of the problem for LSQ approximation of the given data points (x(k); f(k)), x(k) = (x₁(k), x₂(k), x₃(k),), i.e. CFN minimizes the functional:where F(CFN, x) = F(x) is an element of the LSQ approximation.

We define the function F(x₁, x₂, x₃) within the intervals [a_i, b_i], [c_i, d_i], i = 1, …, 3; base functions and coefficients as follows:where

CFN (i₁, i₂, i₃) are elements of coefficient matrix with a size n₁ × n₂ × n₃;

Base functions are .

For the base function we use:

• Algebraic functions: 1, s, s², …, s^k, …;

• Tchebishev’s functions: T₀(s) = 1, T₁(s) = s, …, T_k(s) = 2sT_k-1(s) − T_k-2(s), for k = 2,3, …, i.e. T_k(s) = cos(k arccos(s));

• Trigonometric functions: 1, sin(s), cos(s), …, sin(ks), cos(ks), …;

In our case the base function and coefficient, defined above, are:

References

All Figures

	Fig. 1. One measured Ne profile (red curve), the inferred O+ distribution (light gray line), the transition height (yellow line), and the H+ distribution above transition height (green dashed line). The TSMP H+ distribution is plotted with the dark gray line. The blue curve represents the expected light-ion (H+ + He+) distribution.
In the text

Fig. 2. (a) Daytime Ne profiles measured at three different geomagnetic latitudes: equatorial (green), midlatitude (red), and higher midlatitude (blue). (b, c, and d) The figures in each of the three panels reproduce each of the profiles and the inferred ion distributions in the same format as in Figure 1. The two purple arrows in panel c formulate the difference between the inferred H+ and O+ densities at the transition height as evidence for the existence of extra ionization. Single purple arrows in panels b and d indicate coincidence between the two values.

In the text

	Fig. 3. The difference ln(O+) − ln(H+), for all ISIS-1 profiles, taken at the transition height as a function of geomagnetic latitude is shown on the left side. The geomagnetic latitude/local time contour map of the difference of O+ and H+ densities at the transition height is given at the right panel.
In the text

	Fig. 4. The scatter plot of the modeled g (gm) versus measured g values. The thick blue line is the linear regression through the origin over gm values.
In the text

	Fig. 5. The projection of the function gm on each of the axes plane approximated with low order polynomials; top panel: the linear fit shows that gm have a tendency to decrease with increasing O+ density; middle panel: gm has a minimum of around 0.6 at low latitudes and increases toward higher latitudes; bottom panel: local time dependence shows a minimum in the early afternoon hours.
In the text

Fig. 6. A sample profile measured on day 97 of the year 1969, at 17:03 UT around equator (lat = 12° S, long = 12.5° E). The ion density profiles in logarithmic scale obtained without optimization (denoted as “old”) and those with optimization (denoted by “new”) are color coded, along with the measured Ne profile presented with the black line. The “old” profiles are presented by dashed lines, while the “new” profiles are given by solid lines. The total ion density presented the sum of the “new” profiles is given by the orange line.

In the text

	Fig. 7. The distribution of the RMS errors comparing modeled with measured profiles before and after TaD upgrade for the whole selected ISIS-1 profiles in: (a) boxplot format (left panel) and (b) histogram format (right panel).
In the text

	Fig. 8. Scatter plots of the TEC estimates derived from both ISIS-1 and TaD profiles, separately for the previous (left) and the improved (right) TaD versions and two altitude zones, the plasmaspheric (top) and the topside ionosphere (bottom).
In the text

	Fig. 9. Error distributions between the model and the data-derived TEC estimates for all of the cases considered in the paper.
In the text

	Fig. 10. Examples for TaD-derived profiles based at Malvern site. The ISR EDP are denoted with the cross symbol. Red line indicates the modeled O+ profile and blue line indicates the modeled EDP.
In the text

	Fig. 11. The normalized RMSE distribution of model predictions for O+ in respect to the measured profiles.
In the text

	Fig. 12. The normalized RMSE distribution of modeled total electron density in respect to the measured profiles.
In the text

	Fig. 13. The distribution of the simple difference between the observed ISR and the modeled electron densities at 400 km (top), 500 km (middle), and 600 km (bottom).
In the text

	Fig. 14. The scatter plot of TEC modeled versus TEC observed estimates.
In the text

	Fig. 15. The TEC error (abs(TECobs-TECmod)) distribution and the relative TEC error distribution. From the top to the bottom: (a) results based on all available measurements; (b) the corresponding distribution diagrams for a subset of measurements that correspond to daytime EDP; (c) the distribution diagrams for nighttime cases.
In the text

	Fig. 16. (a) Five Ne profiles in semi-logarithmic scale (ln(Ne), h) measured at low-midlatitudes in afternoon hours. Arrows indicate the corresponding transition heights. (b) The augmented lower part, demonstrating the large variability of the scale height at the topside F region.
In the text

	Fig. 17. A statistical representation of the measured derivatives (measured values marked by red crosses) of the profiles from the whole database. The histogram of Hp obtained from the corresponding profiles is shown at the top portion of the figure.
In the text