Adjustments of the TaD electron density reconstruction model with GNSS-TEC parameters for operational application purposes

Validation results on the latest version of TaDmodel (TaDv2) show realistic reconstruction of the electron density profiles (EDPs)with an average error of 3 TECU, similar to the error obtained from GNSS-TEC calculated paremeters. The work presented here has the aim to further improve the accuracy of the TaD topside reconstruction, adjusting the TEC parameter calculated from TaD model with the TEC parameter calculated by GNSS transmitting RINEX files provided by receivers co-located with the Digisondes. The performance of the new version is tested during a stormperiod demonstrating further improvements in respect to the previous version. Statistical comparison of modeled and observed TEC confirms the validity of the proposed adjustment. A significant benefit of the proposed upgrade is that it facilitates the real-time implementation of TaD. The model needs a reliable measure of the scale height at the peak height, which is supposed to be provided by Digisondes. Oftenly, the automatic scaling software fails to correctly calculate the scale height at the peak,Hm, due to interferences in the receiving signal. Consequently the model estimated topside scale height is wrongly calculated leading to unrealistic results for themodeledEDP.TheproposedTECadjustment forces themodel to correctly reproduce the topside scale height, despite the inaccurate values of Hm. This adjustment is very important for the application of TaD in an operational environment.


Introduction
Knowledge of the state of the upper atmosphere, and in particular its ionized part, is very important in several applications affected by space weather, especially the communications and navigation systems that rely on radio transmission.Although it is forecasted that solar cycle 24 will be weaker than the pervious one (Podladchikova & Van der Linden 2011), solar storms are expected to affect seriously conditions in the Earth's ionosphere and plasmasphere.Therefore an advanced method to nowcast the electron density variations is required for the support of satellite operations.
To better characterize the ionosphere and forecast its disturbances over Europe, a data and model infrastructure platform called the European Digital Upper Atmosphere Server (DIAS) has been established in the National Observatory of Athens by a European consortium formed around eight ionospheric stations, and funded by the European Commission.The DIAS system (http://dias.space.noa.gr)operates since 2006 and the basic products that are delivered are real time and historical ionograms from all DIAS ionospheric stations with the results of the automatic scaling archived in SAO format, frequency plots and maps of the ionosphere over Europe based on the f 0 F2, M(3000)F2, MUF and electron density parameters, as well as long-term and short-term forecasting up to 24 h ahead (Belehaki et al. 2005(Belehaki et al. , 2006a(Belehaki et al. , 2007)).The operational capacity of the basic models implemented in DIAS for the forecast of ionospheric disturbances has been recently validated by Tsagouri (2011).An additional product that is important for spacecraft operations at LEO and MEO heights is the height integral of the electron content 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 the TEC parameters 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.DIAS has the possibility to provide this additional service through the implementation of a topside profiler on the European ionosonde network measurements.The Topside Sounders Model assisted by Digisonde (TaD) can be exploited for this purpose.
TaD is based on the Topside Sounders Model (TSM), which is a set of empirical models 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 (see, e.g., Marinov et al., 2004;Belehaki et al. 2006b;Kutiev et al. 2006;Kutiev & Marinov 2007).To approach the problem of developing a topside/plasmaspheric electron density reconstruction model with increased accuracy, Kutiev et al. (2009a) offered analytical formulas for obtaining the shape of the vertical plasma distribution in the topside ionosphere and plasmasphere based on the TSM parameters.This profiler was named Topside Sounder Model Profiler (TSMP), and it models separately the O + and H + density profiles.To obtain the density distribution, TSMP needs specification of the F layer maximum density (N m F2), its height (h m F2), and its scale height (H m ) at its lower boundary.Digisondes parameters can be used for this purpose, to provide the reconstructed electron density profile from the F layer peak up to GNSS orbits.The new model is the TSMP-assisted Digisonde (TaD) model.
The first edition of the model, TaDv1, was proposed by Kutiev et al. (2009a) who developed a method that applies TSMP to EDP obtained by ground-based Digisonde sounders.Digisonde software (Reinisch & Huang 2001) extends the measured bottomside profile in the topside by assuming an a-Chapman shape with constant scale height H m obtained from the bottomside profile at the F layer peak h m F2.However, it was shown that this is just a rough approximation of the real electron density distribution at the topside (Belehaki et al. 2003).The scale height H m is a neutral scale height, and the Digisonde EDP well above h m F2 decreases upward with a gradient corresponding to double H m .Kutiev et al. (2009a) found that the doubled H m is systematically lower than the TSM scale height H T by a factor of 1.25.This factor was found to vary with 5-10% at the different locations.In order to adjust TSMP to Digisonde measurements, Kutiev et al. (2009a) proposed the multiplication of H m by correction factor k = 2.5 and obtained the corresponding transition height h D through the ratio h D = kH m /R T , with the ratio R T taken from TSM.Having obtained the scale height H D = kH m and transition height h D , the new profiler named TaD (TSMP-assisted Digisonde Profiler) extends the measured bottomside EDP to the topside F region and plasmasphere.It should be noted that while TSMP describes the shape of the topside EDP only, TaD provides EDP from h m F2 up to GNSS orbit heights at Digisonde location.Kutiev et al. (2009b) examined the H + scale height by using ISIS-1 topside sounder profiles which extend up to 3500 km.The H + scale height, denoted as plasmasphere scale height H p , was extracted from ISIS-1 profiles as the largest gradient of the profiles above 2500 km, well above transition height.It was found that the ratio H p /H T is close to 16 at equatorial and low latitudes, as it decreases significantly toward poles.Kutiev et al. (2009b) obtained a simple expression of the ratio depending linearly on squared cosine of geomagnetic latitude.In other words, the H p /H T ratio follows strictly the shape of plasmapause.
TaDv1 was widely verified with all available measured profiles and those extracted from another measurements.Belehaki et al. (2009) summarized comparison results between electron density profiles and TEC parameters extracted from TaDv1 model with (a) CHAMP-derived TEC parameters, (b) CHAMP reconstructed profiles (c) ground-based GNSS-derived TEC parameters and (d) profiles reconstructed from RPI/IMAGE plasmagrams.In all cases TaD follows the general trend of plasmaspheric observations derived from the above datasets.Statistically sufficient comparison with GNSS-derived TEC showed an average model error of 1.92 TECU for Athens and 1.0 TECU for Juliusruh data with standard deviation of 5.2 and 2.9 respectively.
In a recent upgrade (TaDv2) TaD provides calculation of O + , H + , and He + density distributions in transition region between topside F region and plasmasphere, extracted from the analysis of the electron density profiles from ISIS-1, and in addition approximates the plasmaspheric scale height as a function of altitude, latitude, local time, and season using an optimization procedure to achieve best fit with the measured profiles (Belehaki et al. 2012).Comparison of the improved TaD with the previous version gives a reduction of the model error of more than two times.
In this paper, further improvements of TaD reconstruction model (TaDv3) are attempted through adjustments of the integrated TaD EDP with ground-based GNSS-derived TEC.This is a demand because TaD uses statistically based relations for the transition height h T and the plasmasphere scale height H p , which could not represent always accurately current conditions.Furthermore, TaD takes as input ionogram parameters f 0 F2 and h m F2 as well as the calculated by the Digisonde software scale height H m , which as explained above reflects the ''Chapman idealization'' of the F layer and in many cases it cannot represent the topside ionosphere which in reality the F layer is not a real ''Chapman-like''.Even with the adjustment of the H m using a constant correction factor, proposed in the first version of the TaD, deviations from real measurements coud be significant, as shown by Belehaki et al. (2009).In addition in some cases the Digisonde ARTIST software provides unrealistally low values of H m , possibly caused by erroneous automatic scaled ionospheric parameters provided by Digisondes.For all these reasons, adjusting the TaD integral to the TEC parameters derived from ground-based GNSS receivers could result in theoretical N e profiles that are closer to reality, and this is what will be investigated in this paper.

Adjustment of TaD parameters to GNSS-TEC
According to TaD formulation, a key parameter for the reconstruction of the electron density profile (EDP) in topside ionosphere and plasmasphere is the Digisonde-derived scale height H m .Therefore a first priority for the real-time implementation of TaD is the quality check of this parameter, as it is extracted from the SAO files archived in DIAS system.Table 1 shows the statistical results for all Digisonde stations, participating in the DIAS network from 2006 to 2010.The total number of H m values accepted for the analysis is 1,163,564 out of 1,342,067 total SAO files.The first column of Table 1 shows the number of H m values in each row (bin).The corresponding H m value is given in the second column.Third column shows the percentage of the numbers in each bin from the total number of values.The last column shows accumulated percentage from the highest bin to the lowest.According to this quality check the values of H m above 40 km are only 34% of all data.In 66% of cases, When H m is not available or it is not accurately or reliably estimated by the ARTIST software, TaD algorithm does not provide reliable results.To overcome this problem, we propose here a method to adjust the TEC derived from TaD to the TEC derived from GNSS RINEX files to compensate the uncertainties in H m value.This is the third version of TaD (TaDv3).
Figure 1 (Kutiev et al. 2009a).Although it assures a statistical reliability of TaD profiles, it cannot guarantee reliable reconstruction of EDP in all individual Digisonde measurements.One way to improve the accuracy of TaD profiling is to compare the EDP integral (bottomside + topside profiles) with GNSS-derived total electron content (TEC), obtained at the location of the ionosonde.By varying the TaD parameters, we can adjust EDP integral to TEC(GNSS) values.This adjustment cannot guarantee that TaD profile represents the reality, but it is one step ahead in improving the accuracy of the TaD profiling based on average relations.
As will be shown below, the profile shape is most sensitive to the changes of topside scale height H D , respectively, to correction factor k. We prepared a procedure, in which we vary subsequently EDP parameters k, H m , h D , and H p and calculate at each step the difference dTEC between the integral of TaD profile and TEC(GNSS).Bottomside profiles provided by Digisonde software are added to TaD integral.Figure 2 (top panel) shows as a sample the difference dTEC = TEC (TaD) À TEC(GNSS) obtained by varying the correction factor k with a step 0.01.In this sample, TaD is applied to Athens Digisonde profile taken at 12:00 on 4 October 2004.dTEC becomes zero at k = 3.25, denoted as k0.As can be expected, TEC(TaD) is lower than TEC(GNSS) for smaller k values and grows as k increases.It is a common feature of the adjustment procedure and dTEC depends linearly on k around k0. Figure 2 (bottom panel) shows part of the profiles with default value 2.5 (no adjustment -blue curve) and k0 = 3.25 (with adjustmentred curve) for the same case as in the top panel of Figure 2. Note that the change of correction factor k changes also transition height h D and plasmasphere scale height H p .
Adjustment procedure was applied to the hourly Digisonde profiles obtained at Athens and Juliusruh in the 6-year-period 2001-2005.Histograms of the correction factor k are shown in Figure 3.The problem of very low values that the Digisonde software often provides is highlighted here.These values seem not to be realistic.These small values require large correction factor to adjust the corresponding TaD profile to measured TEC(GNSS).In Figure 3  Athens, than at Juliusruh.We have inspected individual ionograms for which ARTIST calculates such low values of H m , and it is obvious that this occurs due to local interference that results in incomplete ionogram traces.The distribution shown in Figure 3 indicates that the most probable value of k is around 2.5, which confirms the previous statistical estimates for the correction factor.TaD EDP is based on seven parameters: N m F2, h m F2, H m , k, R T , h D , and H p .The F layer peak parameters N m F2 and h m F2 are measured, H m is calculated by assuming a-Chapman distribution around the peak, and the rest are statistically derived from the topside sounder data.The statistical relations of k, R T , h D = kH m /R T , and H p = kH m (9 cos 2 (glat) + 4), obtained by Kutiev et al. (2009b), introduce errors (uncertainties) even when TEC(TaD) is adjusted to TEC(GNSS) with a good accuracy.To estimate the errors coming from all parameters in the TaD profile during adjustments, we take the profile calculated with correction factor k0 and subsequently (one by one) start varying k, H m , h D , and H p around their values at k0.We consider the total change DT of TaD integral, resulting from variation of profile parameters, as a total derivative composed by partial derivatives of T on individual parameters at the zero value k0, multiplied by their respective shifts: The partial derivatives are obtained by varying parameter values around their initial values at k0 and numerically calculate respective derivatives at the initial values.If we take the shifts Dk, DH m , Dh D , and DH p to be 10% of their values at k0, the terms in Eq. ( 1) have the following contributions in percent: The first term of Eq. ( 1) contributes with more than 90% to the total change of integral T. When all parameters of the profile change at 10% from their initial values, the largest part of total change of T is due to k change.It means that the integral of TaD EDP is most sensitive to the changes of correction factor k. If we regard H m as an external parameter that virtually does not contain error, k remains the only parameter controlling the integral of TaD profile.From this result we can derive an important advantage of the proposed adjustment procedure: we can make adjustment by solely varying correction factor k.
Figure 4 shows oT =ok versus TEC (TaD) for all data from Athens (bottom) and Juliusruh (top).Both dependences are pretty linear and the partial derivative can be easily obtained from the linear regressions given in the plots.If TaD EDP is routinely calculated with default k = 2.5 and then its integral compared with TEC(GNSS), the difference dTEC can be written as: and from here: where TEC 2.5 is TEC(TaD) with k = 2.5, TEC GNSS is the measured TEC(GNSS), and oT =ok is calculated from TEC(TaD)  from regression equation for the considered station.When the proper correction factor k0 is obtained, adjusted TaD profile can be obtained with one more run.This simple procedure can avoid numerous runs with varying k, iteratively searching dTEC = 0.It should be noted that TaDv3 does not have any threshold for deciding whether the scaling (or correction) factor k needs to be applied.The correction factor is applied to all data.

The performance of the TaD adjusted algorithm
Because Digisonde H m is frequently very low (see Table 1), in addition to the TEC adjustment, we propose the following limits to the Digisonde data in order to ensure the reliable operation of TaD: d transition height h T not less than h m F2 + 100 km; this threshold in the difference between h m F2 and h T is more an empirical conclusion rather than a statistical result, d the difference TEC(GNSS) -TEC(TaD) to go through zero when varying k, or at least the closest approach to be less than 0.3 TECU.
A test run of TaD with real-time SAO files was done during a storm period.The following combinations of the model versions were tested: d TaDv3t: this is the third version of TaD with no limits applied to Digisonde data, d TaDv3e: this is the third version of TaD with the abovedefined limits applied to Digisonde data, d TaDv1e: this is the first version of TaD with the abovedefined limits applied to Digisonde data.
It is important to clarify here that the TaDv3 includes the optimization procedure adopted in TaDv2 (see Belehaki et al. 2012).
The calculated TEC values with all above combinations, during October 20-26, 2001 storm, are presented in Figure 5, together with the TEC(GNSS) values over Athens and Juliusruh.For the application of TaDv1 a constant correction factor k = 2.5 was used, while for the application of TaDv3 the k factor was adjusted according to the technique described in Section 2. As one can see TaDv3 (both options) performs much better than TaDv1.
To further investigate the impact of low h m F2 values that Digisonde often provides, we present in Figure 6 the scatter plots between TEC(GNSS) and TEC(TaD) over Athens and Juliusruh.Here we run the TaDv3 with the above limits applied (t4, red crosses) and with no limits applied (e4, purple diamonds).When no limits are applied, there is a large scatter between TEC(TaD) and TEC(GNSS), with TEC(TaD) < TEC(GNSS).This means that TEC(TaD) cannot be adjusted to reach the TEC(GNSS) values, when k is varying between 1.5 and 6.5.The green line is a regression fit to e4 values.
This shows clearly that the above limits are valid and necessary to be applied for operational purposes.

Discussion and conclusions
The paper describes a method for improving the accuracy of TaD EDP by adjusting its integral to the TEC derived from GNSS signals (TaDv3).The aim of the paper is to describe the theoretical backgrounds of TaDv3, to assess the possible errors, and to demonstrate its performance.TaD EDP of the topside and plasmasphere accomplish the bottomside profile provided by Digisonde software.The EDP above the F layer peak is composed by the O + profile, representing the topside F region and H + profile representing the plasmasphere, both equalized at the upper transition height.In TSMP, the transition height h T is connected through the ratio R T to the topside scale height H T (h T = H T /R T ); the plasmasphere scale height H p is also expressed by H T and depends on geomagnetic latitude.Therefore, the whole EDP above the peak is defined by H T and R T through statistical relations provided by TSM.TaD couples TSMP to Digisonde profiling technique, by substituting H T with Didisonde scale height H m and the correction factor k. In TaDv3 profiler, the key parameters become the correction factor k and the TSM ratio R T .It is assumed that the ratio of topside scale height and transition height is an inherent parameter of topside ionosphere and when the topside scale height is changed, the corresponding transition height is changed accordingly.The ratio R T is not a constant throughout the globe; TSM provides it as a function on time of the year, local time, geomagnetic latitude, solar flux F 10.7 , and K p .Therefore, for each individual calculation of EDP, TaD is provided by respective value of R T .
It is shown that the correction factor k varies significantly when TaD EDP integral is adjusted to TEC(GNSS).Adjusting of TaD EDP integral (TEC(TaD)) to the measured TEC (TEC(GNSS)) is accomplished by varying k between 1.5 and 6.0.In this range, the TEC difference (dTEC) is pretty linear with k, which simplifies the iteration procedure.The value of k for dTEC = 0, denoted as k0, is statistically obtained over 6 years of data from Athens and Juliusruh.The factor k0 compensates the lower values of H m in respect to the topside scale height and its most probable value is around 2.5, as it was found in the previous studies.
As was pointed out above, adjustment of integral of TaD EDP to the measured TEC does not guarantee that the profile itself is unique.The weaker statistical dependences, mentioned above, introduce errors, which we tried to estimate.The possible errors coming from each parameter of the profile are estimated by the change of the integral when the parameter value is varied in certain range.To estimate the total error when varying individual parameters, we consider the total change T of the integral as a total derivative DT composed of partial derivatives of T on individual parameters multiplied by their respective shifts.It was found that changes of correction factor k contribute to more than 90% of the total change of the integral.This is an important clue for further simplification of adjusting procedure: TEC(TaD) can be changed by solely changing the correction factor k.
To check the validity of the proposed adjustment and the improvement of EDP accuracy, we need to have measured EDP in the topside F region and plasmasphere.Unfortunately, the topside sounders operated well before Digisonde groundbased sounders and direct comparison is not possible.One way for partial verification of the proposed procedure is to compile simultaneous measurements at one location of topside EDP from incoherent scatter radar (ISR), Digisonde EDP and scale height, and GNSS-derived TEC.Such verification will be a subject for a follow-on paper.Fig. 6.The scatter plots between TEC(GNSS) and TEC(TaD) over Athens (at the top) and Juliusruh (at the bottom) for the two versions of the optimized TaD: e4 (purple diamonds) and t4 (red crosses).The green line is a regression fit to e4 values.
illustrates the difference between Digisonde and TaD EDPs.Digisonde profile (red curve) is obtained from Ebro ionosonde station at 07:00 UT on 5 June 2001.Digisonde scale height H m = 49.8 km, N m F2 = 1.0 • 10 6 cm À3 , and h m F2 = 284.3km.The figure shows O + and H + profiles, which compose TaD EDP.TaD topside scale height H D = 124.5 km, H p = 1490 km, and transition height h D = 825 km.The correction factor k, which converts H m into H T , is statistically obtained by comparing Digisonde data with the corresponding TSM product

Fig. 1 .
Fig. 1.The difference between Digisonde and TaD EDPs.Digisonde profile (red curve) is obtained from Ebro ionosonde station at 07:00 UT on 5 June 2001.The figure shows O + (magenta line) and H + profiles (dashed blue line), whose sum compose TaD EDP (green line).The transition height is indicated with the yellow line.

Fig. 2 .
Fig. 2. Top panel: the difference dTEC = TEC(TaD) -TEC(GNSS) obtained by varying the correction factor k with a step 0.01.In this sample, where TaD is applied to a profile from Athens Digisonde, dTEC becomes zero at k = 3.25, denoted as k0.Bottom panel: the same profile case shown in the top panel, shown here with default value k = 2.5 (no adjustment -blue curve) and k0 = 3.25 (with adjustment -red curve).

Fig. 3 .
Fig. 3. Adjustment procedure was applied to the hourly Digisonde profiles obtained at Athens and Juliusruh in the 6 years period 2001-2005.Histograms of the correction factor k are given separately for all H m values (red bars) and for those greater than 40 km (open blue bars).

Fig. 5 .
Fig. 5.The TEC GNSS and TEC TaD values during October 20-26, 2001 storm over Athens (at the top) and Juliusruh (at the bottom).TEC TaD values were obtained for both TaD versions: TaDv1 with constant correction factor k = 2.5 and TaDv3 with adjusted correction factor.

Table 1 .
Number of Digisonde-derived H m values in each 20 kmwide bins in altitude, persentage from the total number and accumulated percentage in each bin.The numbers in the first column show the lower value of the bin range.
H m drops to very low values (<40 km), which is unrealistic.In statistical analysis, we take the value of 40 km as a boundary between the higher H m which perfectly cooperate with TaD and those below which may lead to unrealistic profiles.The value of 40 km was chosen because its doubled value (80 km) is approximately one standard deviation (40 km) below the average value of TSM H T (120 km).H m higher than 40 km would yield topside scale height H D within the normal distribution of H T .This boundary is considered only for statistical purposes, to demonstrate how H m influences the TaD performance.