Three-dimensional data assimilation for ionospheric reference scenarios

The reliable estimation of ionospheric refraction effects is an important topic in the GNSS (Global Navigation Satellite Systems) positioning and navigation domain, especially in safety-of-life applications. This paper describes a three-dimensional ionosphere reconstruction approach that combines three data sources with an ionospheric background model: spaceand ground-based total electron content (TEC) measurements and ionosonde observations. First the background model is adjusted by F2 layer characteristics, obtained from space-based ionospheric radio occultation (IRO) profiles and ionosonde data, and secondly the final electron density distribution is estimated by an algebraic reconstruction technique. The method described is validated by TEC measurements of independent ground-based GNSS stations, space-based TEC from the Jason 1 and 2 satellites, and ionosonde observations. A significant improvement is achieved by the data assimilation, with a decrease in the residual errors by up to 98 % compared to the initial guess of the background. Furthermore, the results underpin the capability of space-based measurements to overcome data gaps in reconstruction areas where less GNSS ground-station infrastructure exists.


Introduction
The ionosphere is the upper part of the atmosphere where sufficient free electrons exist to affect the propagation of radio waves.Therefore, the ionospheric parameters, such as three-dimensional electron density distribution, the iono-spheric layer peak characteristics, and the total electron content (TEC), are important information for Global Navigation Satellite Systems (GNSS) users.The ingestion of actual ionospheric measurements into a background model, such as NeQuick (see Nava et al., 2008) or International Reference Ionosphere (IRI) (see Bilitza, 2001;Bilitza and Reinisch, 2008), is a commonly applied technique for estimating the ionospheric parameters.Several approaches had been tested for the ionospheric imaging combining actual direct and indirect measurements with either an empirical or a physical model background (see e.g.Angling and Cannon, 2004;Angling, 2008;Schunk et al., 2004;Wang et al., 2004;Scherliess et al., 2009;Brunini et al., 2011;Pezzopane et al., 2011;Galkin et al., 2012;Minkwitz et al., 2015Minkwitz et al., , 2016;;Gerzen and Minkwitz, 2016).
Within this work we present a two-step three-dimensional reconstruction approach, which firstly adjusts the initial background model NeQuick (see Nava et al., 2008) by the assimilation of F2 layer characteristics, obtained from space-based ionospheric radio occultation (IRO) profiles and ionosonde data.Secondly, the final electron density distribution is estimated by the algebraic reconstruction technique SMART+ (Gerzen and Minkwitz, 2016) assimilating spaceand ground-based TEC measurements as well as ionosonde observations.The described method is validated by slant TEC measurements from independent ground-based GNSS stations, vertical TEC from the Jason 1 and 2 satellite and ionosonde observations.Special attention is paid to one question: How strong is the potential improvement achievable by the preconditioning of the background with space-based and ionosonde measurements compared to the usage of groundbased TEC only?The approach presented is applied in the scope of the project DAIS (Data Assimilation Techniques for Ionospheric Reference Scenarios) to generate synthetic ionospheric reference scenarios (IRSs) for the validation of the European Geostationary Navigation Overlay Service (EGNOS).EGNOS provides value-added services for the estimation of ionospheric refraction effects.The IRSs are introduced by the ESA in order to conduct the EGNOS end-to-end performance simulations and to assure the integrity of the EGNOS system and associated services (see Arbesser-Rastburg, 2004;Schlüter et al., 2013).
The paper is organized as follows.At first the chosen reconstruction area and periods as well as the applied database are described.Section 3 then explains the reconstruction approach.Section 4 presents the validation results, and finally, the results are summarized and discussed.
2 Database and data processing

Reconstruction area and validation periods
The reconstruction approach described is tested over the extended EGNOS V3 service region (−100 to 110 • E in longitude and −90 to 90 • N in latitude) with a 5 min cadence.We select two periods in 2011, a calm period between the days of year (DOYs) 009 and 026 (9 to 26 January) and a storm period between DOY 286 and 303 (13 to 30 October).The selection of these days is based on the solar radio flux index (F10.7),the two geomagnetic indices Kp and Dst, and the along-arc total electron content rate AATR RMS1hour (see Schlüter et al., 2013;Juan et al., 2014) The impact of such conditions on the EGNOS system is given in Fig. 2. It shows the availability of the EGNOS service for aviation approaches with vertical guidance (APV-1) on DOYs 010 (left) and 298 (right) of the year 2011 (see https://egnos-user-support.essp-sas.eu/new_egnos_ops/apv1_availability) where the EGNOS service availability on DOY 298 is crucially limited.

Ground-based GPS TEC measurements
We use absolute TEC measurements from ground-based IGS stations as input for the assimilation and validation.The details of the absolute TEC calculation are given in Jakowski et al. (2011).For this study, the 1 Hz GPS data of the IGS   1, are excluded from the reconstruction procedure and considered as independent references for the validation.

Ionosonde measurements
The F2 layer characteristics, in particular the critical frequency, foF2, and the peak height, hmF2, of 48 ionosonde stations are collected from the SPIDR database.Their corresponding locations are depicted in Fig. 3.For validation purposes, the data of seven ionosonde stations are excluded from the reconstruction.Furthermore, the measurements of ionosonde stations JI91J and KJ609 are used both for the assimilation and the validation.All stations used for validation purposes are listed in Table 2.
In the following, the NmF2 values, calculated from the foF2 measurements by NmF2 [m −3 ] = 10 10 • 1.24 • foF2 [MHz] 2 , are also denoted as ionosonde measurements.In order to avoid artefacts in the data assimilation procedure, the ionosonde data are filtered before the application.For more details the reader is referred to Gerzen et al. (2015).

Space-based GPS TEC measurements and ionospheric radio occultation profiles
For both periods, the data of the low Earth orbit (

TEC derived from altimeter data
Dual-frequency altimeter missions, such as Jason 1 and 2, are an excellent source for vTEC data independent of the GNSS systems.In contrast to the GNSS-derived vTEC, the altimeter measurements are naturally vertical.

Scenario reconstruction approach
Figure 4 outlines a sketch of the developed assimilation process.The single steps are further detailed in the subsections of this section.
The correct characterization of the vertical shape of the profiles becomes a difficult task when assimilating only ground-based sTEC because of limited vertical information included in these data (see e.g.McNamara et al., 2008McNamara et al., , 2011;;Minkwitz et al., 2015;Gerzen and Minkwitz, 2016).The inclusion of space-based sTEC improves the geometry situation.However, the adjustment of the background in terms of the F2 layer characteristics before starting the assimilation procedure seems to be especially advisable (e.g.Bidaine and Warnant, 2010) since the F2 layer dominates the shape of the whole profile.
Thus, we first estimate global NmF2 and hmF2 maps.For that purpose, the modified successive correction method (MSCM; see Gerzen et al., 2015) is applied.By means of MSCM, F2 layer characteristics from ionosondes and IRO profiles are assimilated into the corresponding twodimensional background models.As global background models, the Neustrelitz Peak Density Model (NPDM; see Hoque and Jakowski, 2011) and the Neustrelitz Peak Height Model (NPHM; see Hoque and Jakowski, 2012) are deployed.Subsequently, the NeQuick model is adjusted by the incorporation of these maps.
The adjusted version of the NeQuick model serves as the initial guess for SMART+ (see Gerzen and Minkwitz, 2016), which is a fusion of the algebraic iterative tomography method SMART and a three-dimensional successive correction method (SCM).SMART+ assimilates the groundand space-based TEC as well as F2 layer characteristics from ionosondes and IRO profiles available over the observed area.SMART distributes the available integral observations to the electron densities of the voxels intersected by at least one ray path, whereas the three-dimensional SCM introduces spatial correlation between regions covered by measurements and those not covered.
To give a special focus on the question of how the inclusion of additional measurements (additional to the groundbased TEC) and the preconditioning of the background influence the assimilation results, we compare two three-

Two-dimensional assimilation by MSCM
For the assimilation of the F2 layer characteristics, MSCM is used with the same configurations as detailed in Gerzen et al. (2015).Contrary to the version used within Gerzen et al. (2015) though, both the ionosonde and IRO profile data are assimilated by MSCM here.The first estimate of the unknown parameters hmF2 or NmF2 is given by the models NPHM or NPDM respectively.Then MSCM iteratively fits the first estimate towards the measurements in the neighbourhood.
Figure 5 presents as an example NmF2/foF2 (left) and hmF2 (right) maps reconstructed for DOY 295 in 2011 at 12:30 UT.The locations of assimilated measurements are marked magenta (dots for ionosondes and stars for IRO).The global maps are reconstructed with a 15 min cadence and are interpolated to a 5 min cadence by the linear interpolation described in Schaer et al. (1998).

Adjustment of the NeQuick model by the reconstructed F2 layer maps
To include the updated F2 layer characteristics, the internal foF2 and hmF2 calculations of the NeQuick version 2.0.2 are replaced by the reconstructed foF2/NmF2 and hmF2 maps.Moreover, the internal NeQuick model function of M(3000)F2 is replaced by the following simplified relation: M (3000) F2 = 1490 hmF2+176 (see Shimazaki 1955).In this way the whole shape of the NeQuick electron density profiles is adjusted.
Figure 6 illustrates the possible improvement of the adjustment described.Electron density profiles are calculated at the position of the independent ionosonde station TO536 (see Table 2) on DOY 295.The NeQuick model profile is depicted in green; the profile calculated after the inclusion of the reconstructed maps, only assimilating ionosonde observations, in blue; the adjustment with ionosonde and IRO observations in red; and the measurement of the ionosonde TO536 as a violet dot.

Three-dimensional assimilation by SMART+
The TEC is related to the electron density Ne by TEC = Ne (h, λ, ϕ) ds, where s is the corresponding ray path between GNSS satellite and receiver and Ne (h, λ, ϕ) is the unknown electron density function depending on altitude h, geographic longitude λ, and latitude ϕ.
We discretize the ionosphere by a voxelized threedimensional grid with a horizontal spatial resolution of 2.5 • .The altitude resolution is 30 km for altitudes between 60 and 1000 km and increases exponentially with increasing altitude for altitudes above 1000 km.In total, 54 altitude steps are obtained.Assuming the electron density function to be constant within a voxel, we derive a linear system of equations (LSE): where y is the vector of the TEC measurements, x i = Ne i is the electron density in the voxel i, and a si is the length of the ray path s in the voxel i.
The SMART+ method (see Gerzen and Minkwitz, 2016) is applied here to solve the LSE.SMART+ combines SMART and three-dimensional SCM.First, the iterative algebraic tomography method SMART starts with the initial guess provided either by the NeQuick model (version A) or the adjusted NeQuick model (version B).Within the iteration step of SMART, the multiplicative innovation is given by a weighted mean of the deviations between the measurements and the current estimate of the measurements.The electron densities of the voxels not intersected by at least one ray path remain equal to the initial guess.
In order to assimilate the F2 layer characteristics using SMART, NmF2 is expressed as TEC: for each measured NmF2-hmF2 pair the corresponding voxel number q is calculated from the corresponding latitude, longitude, and hmF2 information.Then the corresponding TEC value is defined as TEC = NmF2 • s q , where s q is the length of the voxel q in the longitudinal direction.
The iteration process stops either after a predefined number of iteration steps (we used 26) or if the mean TEC deviation from the assimilated measurements goes below a predefined threshold.Thereafter, an extrapolation from the intersected voxels to those not intersected by any TEC ray path is done using the three-dimensional SCM method, assuming a Gaussian covariance model for the electron densities (see Gerzen and Minkwitz, 2016).Figure 7 presents as an example the reconstructed electron density on DOY 295, 12:30 UT.

Calculation of the IRSs
The IRSs, i.e. vTEC maps, are calculated from the three-dimensional reconstructions by an integration of the electron density profile values using vTEC(λϕ) = 53 i=1 Ne rec (altitude i , λ, ϕ) • (altitude i+1 − altitude i ). Figure 8 presents an example IRS calculated from the reconstructed electron density shown in Fig. 7.

Validation of the adjusted NeQuick with ionosonde characteristics
For the validation of the adjusted NeQuick (Step 1), the NmF2 and hmF2 data of the ionosondes listed in Table 2 serve as references.At these stations, the electron density profiles are calculated with 10 km altitude resolution.
Here, three versions of the electron density profile output are compared: one is calculated with the NeQuick model, the second is calculated with the NeQuick version assimilating only ionosonde observations in Step 1, and the third is the adjusted NeQuick assimilating ionosonde and IRO data.Afterwards, we derive from these profiles (separately for each of the three versions) the F2 characteristics, i.e.NmF2 reconstructed and hmF2 reconstructed .Then we derive the residuals: NmF2 measured − NmF2 reconstructed and hmF2 measured − hmF2 reconstructed are calculated.Thereby NmF2 measured and hmF2 measured denote the characteristics measured at the reference ionosonde stations (see Table 2).In Fig. 9, the distribution of the residuals over all reference ionosonde stations is shown, including the mean, standard deviation (SD), and RMS of the residuals.Again, the colour green is used for the results calculated with the NeQuick model, blue for the NeQuick model assimilating solely ionosonde data, and red for the adjusted NeQuick model assimilating ionosonde and IRO data.
The residual statistics indicate the advantage of the preconditioning of the background using the current F2 layer characteristics.The mean of the hmF2 residuals is decreased by up to 98 % and the RMS value by up to 27 % compared to the corresponding pure NeQuick residual statistics.The NmF2 statistics are also up to 81 % lower (in relation to the magnitude).The inclusion of the IRO data in the assimilation procedure induces a slight additional improvement during the storm period.For the quiet period, no additional improvements are visible.Therefore, one of the reasons might be inconsistencies between the ionosonde and IRO data.During the quiet period, the deviations between the hmF2 and NmF2 values calculated by the NeQuick model and the measurements are comparatively small.Thus, even small inconsistencies in the data can significantly influence the assimilation and validation results.

Validation of the IRSs with altimeter data
In this subsection we compare the absolute residuals between the Jason 1 and 2 vTEC measurements and the reconstructed vTEC derived from the reconstruction versions A and B, in detail: |vTEC jason − vTEC ground | and |vTEC jason − vTEC all | respectively.Since the Jason satellites fly at an orbit height of about 1336 km, only voxels up to this height are integrated to calculate vTEC ground and vTEC all .In Fig. 10, the distribution of the absolute residuals is given for the two investigated periods in five TECU (total electron content unit) bins (the last bin sums up all the higher values).The vTEC ground residuals are depicted in green, the vTEC all in red.The majority of all residuals have values less than 10 TECU.
Table 3 presents the mean, SD, and RMS values of the absolute residuals of Fig. 10.Version B decreases the mean values up to 0.4 TECU compared to version A. However, in most cases the SD of the version B residuals is higher than for version A. We assume that the partial increase in the residual deviation may be caused by inconsistencies in the assimilated ground-and space-based sTEC (see Sect. 4.3).Furthermore, there may be inconsistencies between the GPS and altimeter vTEC (see Azpilicueta and Brunini, 2009).
The histograms and statistics of the vTEC ground and the vTEC all residuals are similar.Thus, to clarify the differences between versions A and B, the ratio |vTEC jason −vTEC all | |vTEC jason −vTEC ground | is additionally considered.This quotient is smaller than one if the additional inclusion of space-based and ionosonde data introduces an improvement.Figure 11 (quiet) and Fig. 12 (storm) present the percentage proportion between the numbers of quotients that are smaller than one and the total number of samples on 1 day.In other words, the number of improvements in percent due to inclusion of space-based and ionosonde data is shown.The green line indicates the 50 % value.The number of improved residuals is around 60 % for most of the days within the investigated periods (∼ one-third higher than 60 %).The improvements are more visible during the storm period.

Validation of the three-dimensional reconstructions with ground-based sTEC
To assess the capability for estimating sTEC, the following parameters are compared: sTEC model derived from the initial three-dimensional electron densities calculated by the NeQuick model, sTEC ground derived from the reconstruction version A, and sTEC all from reconstruction version B are separately compared to the measurements of the reference  stations (see Table 1).For each period the residuals between the reconstructed and the measured sTEC values are calculated as dTEC all = sTEC measured − sTEC all , dTEC ground = sTEC measured −sTEC ground , and dTEC model = sTEC measured − sTEC model .
Figure 13 exemplarily depicts the histograms of the sTEC residuals at the station Kiruna.On the left-hand side, the distribution of the residuals for the quiet period is shown and on the right-hand side the disturbed period is shown.An overestimation by the NeQuick model (green) is visible for both periods in accordance with the results presented in Nigussie  The measured sTEC ranges up to 140 TECU during the quiet period and up to 320 TECU during the storm period.For the quiet period, improvements from sTEC model to sTEC ground and from sTEC ground to sTEC all are observable and are underpinned by the correlation coefficients.The corresponding P value is around zero, indicating a high significance for the correlation coefficient.During the disturbed period, the behaviour of sTEC model and sTEC ground is very similar.Conversely, for sTEC all we observe that several smallmagnitude outliers of the NeQuick model (see Fig. 14, right column, third row) are removed.However, several highmagnitude outliers are simultaneously introduced (Fig. 14, right column, first row).These outliers are not visible in the sTEC ground figure (Fig. 14, right column, second row).Thus, they are most probably introduced by inconsistencies between the assimilated space-based data and the reference ground-based sTEC.

Validation with the ionosonde data
The validation of the reconstructed electron density profiles with the reference NmF2 and hmF2 data shows a significant decrease in the evaluated statistics compared with the pure NeQuick results.The decrease is up to 98 % for the hmF2 residuals and up to 81 % for the NmF2 residuals.The inclusion of the IRO profile parameters in the preconditioning procedure causes an additional improvement during the storm period.

Validation with Jason data
Both versions A (assimilation of ground-based data only) and B (assimilation of ground-, space-based, and ionosonde data) are validated with the Jason 1 and 2 vTEC measurements.During the storm period, this comparison indicates a slightly lower mean value of the absolute residuals of version B compared to version A.

Validation against ground-based sTEC
The assimilation of ground-based sTEC clearly improves the initial guess of the NeQuick model.The median of the TEC residuals is decreased by up to 65 % during the quiet period and by up to 62 % during the storm period.However, no or only a very small advance (up to 16 %) is observed after the additional inclusion of the space-based and ionosonde data in the assimilation procedure.This is probably due to inconsistencies between the assimilated space-based data and the reference ground-based sTEC as detailed in the validation section.

Conclusions and discussion
Within this work, time series of three-dimensional electron density and IRSs, representing quiet and disturbed ionospheric conditions, are generated and cross-validated.These electron density values are deduced from the threedimensional assimilation of ground-and space-based TEC and F2 layer characteristics into the NeQuick model.
The validation results show that a crucial improvement is achieved by the adjustment of the whole background to the measured F2 layer characteristics within a preconditioning step.A decrease in the hmF2 residual statistics of up to 98 % is hereby observed.For this preconditioning step, the IRO data are especially important due to the limited availability of ionosonde measurements, in particular over the ocean and in Africa.
Through the subsequently assimilation of the ground-and space-based TEC into the preconditioned background a decrease of the sTEC residual statistics up to 62 % (Kiruna station) is gained.The space-based sTEC measurements cover wide regions where ground-based data are sparse (see Fig. 15).Hence, the use of additional LEO satellites (like SWARM and GRACE) looks very promising to fill the remaining data gaps.A further advantage of the space-based data is their measurement geometry, which is complementary to the angle-limited geometry of the ground-based TEC data.
Summarizing the results, on the one hand, the validation clearly shows the potential of data assimilation especially when combining the data from different data sources,.On the other hand, several results of this study indicate inconsistencies in the data obtained from different measurement instruments.This reflects a serious problem for all data-driven approaches and dealing with such inconsistencies is still a challenging task.Advanced filter methods, allowing the simultaneous filtering of all used data types, may contribute significantly to the solution for this problem.
Moreover, the inverse problem behind the ionosphere reconstruction remains ill-posed and highly underdetermined with angle-limited measurement geometry, even if all currently available measurements are used.Due to the sparseness of the available data and the limited geometry, potential artefacts in the final electron density reconstructions cannot be ruled out and have to be treated carefully.This underlines that the appropriate estimation of the correlation, between the electron density at a measurement location and a location where no measurements are available, is a necessity for an improved reconstruction of the ionosphere.We are currently working on methods that estimate the correlation lengths and variance components of parametric covariance models of the electron densities (see Minkwitz et al., 2015Minkwitz et al., , 2016)).However, a systematic study that investigates the spatial and temporal electron density correlations is highly recommended.

Figure 1 .
Figure 1.AATR box plot for the nya2 station.Values for the disturbed period are presented in the top panel, while calm period values are in the bottom panel.
. The AATR RMS1hour values are calculated at each GNSS ground-based station of the International GNSS Service used (IGS, see Sect.2.2).The Kp and Dst indices are downloaded from the Space Physics Interactive Data Resource (SPIDR) of NOAA's National Geophysical Data Center and the World Data Center for Geomagnetism (WDC) Kyoto.During the calm period, F10.7 was in the range of 75 to 90 solar flux units, Kp index was below 4, and Dst was between −40 and 40 nT.Conversely, the storm period reveals increased solar and geomagnetic activity, with a F10.7 of 120-170 flux units and a severe geomagnetic storm observed during the DOYs 297-298 with a Kp index above 7 and Dst values below −130 nT. Figure 1 presents as an example the AATR RMS1hour response to the ionospheric conditions at the IGS ground-based station nya2 in Norway.The AATR RMS1hour values of each day are organized as a standard box plot.As expected, the top panel shows higher values during the disturbed period than during the calm period depicted in the bottom panel.

Figure 3 .
Figure 3. Ionosonde station positions.Black dots are ionosonde stations used only for assimilation, red triangles are ionosonde stations used only for validation, and black triangles inside red triangles are ionosonde stations used for assimilation and validation.

Figure 6 .
Figure 6.Electron density profiles: NeQuick model (green) and adjusted NeQuick model (blue is assimilation of ionosondes only, red is assimilation of ionosonde and IRO F2 layer data) compared with the independent ionosonde measurement (violet dot).

Figure 8 .
Figure 8. IRS calculated from the reconstructed three-dimensional electron density.

Figure 10 .
Figure 10.Histogram of the vTEC absolute residuals for Jason 1 (left) and 2 (right) satellites during the quiet (top) and disturbed (bottom) periods.

Figure 11 .
Figure 11.The percentage number of reduced absolute Jason residuals per day of the quiet DAIS period.

Figure 12 .
Figure 12.The percentage number of reduced absolute Jason residuals per day of the disturbed DAIS period.

Figure 14 .
Figure 14.sTEC measured versus sTEC all (top row), sTEC ground (middle row), and sTEC model (bottom row) at the cro1 reference GNSS station in the USA during the quiet (left column) and disturbed (right column) periods.

Figure 15 .
Figure 15.Coverage of the reconstructed array by space-based data (right) and by ground-based data (left) exemplarily on DOY 011 in 2011, 09:30 UT.A pixel is coloured blue if at least one voxel above this pixel is intersected by at least one sTEC ray path.

Table 1 .
IGS stations that provide independent sTEC measurements for the validation.

Table 2 .
Ionosonde stations used for validation purposes.Ionosonde stations with italic font are used for assimilation and validation.Station ID Country Lat. ( • N) Long.( • E)

Table 3 .
Mean, RMS, and SD of the absolute Jason residuals.

Table 4 .
The statistics of the absolute sTEC residuals in TECU for DOYs 009-021.

Table 5 .
The statistics of the absolute sTEC residuals in TECU for DOYs 286-303.