Ionospheric tomography by gradient-enhanced kriging with STEC measurements and ionosonde characteristics

The estimation of the ionospheric electron density by kriging is based on the optimization of a parametric measurement covariance model. First, the extension of kriging with slant total electron content (STEC) measurements based on a spatial covariance to kriging with a spatial– temporal covariance model, assimilating STEC data of a sliding window, is presented. Secondly, a novel tomography approach by gradient-enhanced kriging (GEK) is developed. Beyond the ingestion of STEC measurements, GEK assimilates ionosonde characteristics, providing peak electron density measurements as well as gradient information. Both approaches deploy the 3-D electron density model NeQuick as a priori information and estimate the covariance parameter vector within a maximum likelihood estimation for the dedicated tomography time stamp. The methods are validated in the European region for two periods covering quiet and active ionospheric conditions. The kriging with spatial and spatial–temporal covariance model is analysed regarding its capability to reproduce STEC, differential STEC and foF2. Therefore, the estimates are compared to the NeQuick model results, the 2-D TEC maps of the International GNSS Service and the DLR’s Ionospheric Monitoring and Prediction Center, and in the case of foF2 to two independent ionosonde stations. Moreover, simulated STEC and ionosonde measurements are used to investigate the electron density profiles estimated by the GEK in comparison to a kriging with STEC only. The results indicate a crucial improvement in the initial guess by the developed methods and point out the potential compensation for a bias in the peak height hmF2 by means of GEK.


Introduction
The spatial-temporal state of the ionospheric electron density distribution, including, for example, the peak characteristics of the different ionospheric layers and the slant total electron content (STEC) along a ray path, is useful information for almost all radio systems (see Bust and Mitchell, 2008).Tomography is a commonly applied tool to study and monitor these ionospheric key parameters.In this scope, various approaches have been developed.One common idea is the stabilization of the ill-posed inverse problem of ionospheric tomography by means of a physical or empirical background model providing a first guess of the electron density distribution.Subsequently the given electron densities are modified according to the measurements without touching the model coefficients itself (see, e.g., Bust et al., 2004;Wang et al., 2004;Wen et al., 2007;Hobiger et al., 2008;Scherliess et al., 2009;Angling and Jackson-Booth, 2011;Norberg et al., 2015;Gerzen and Minkwitz, 2016).
Since the ionospheric measurements are sparsely available compared to the desired spatial-temporal resolution of the ionospheric state, there is demand to spread out the information to locations adjacent to the measurement geometry.For that purpose, the definition of spatial and temporal correlations between electron densities at different locations is essential and incorrect specifications degrade the estimates of the ionospheric parameters.Motivated by the already successful application of kriging for the estimation of TEC grids/errors for the different satellite-based augmentation systems (Sparks et al., 2011;Sunda et al., 2013) and the provision of regional and global TEC maps (Orús-Pérez, 2005;Sayin et al., 2008), we advance the 2-D kriging to a 3-D kriging approach introduced in Minkwitz et al. (2015).This approach has its origin in the geostatistics (see Wackernagel, 2003;Chilès and Delfiner, 2012), estimates the optimal covariance between two electron densities by the STEC measurements and predicts the electron density based on that.
In general the kriging approach is capable to assimilate various direct and indirect measurements of the electron density, such as ground-based STEC, ionosonde measurements and satellite-based STEC and occultation measurements.In this paper the kriging with STEC measurements based on a spatial covariance model is developed to a kriging based on a spatial-temporal covariance assimilating the STEC measurements of a sliding window.Furthermore, to mitigate the well-known bias in the estimation of the F2 layer peak height hmF2 by ground-based STEC measurements only, the assimilation of ionosonde measurements by means of kriging with derivative information is investigated.In the literature, this method is often referred to as gradient-enhanced kriging (GEK) (see Han et al., 2013;Ulaganathan et al., 2015).The three methods presented start with an initial guess of the electron density provided by an arbitrary 3-D electron density model.In this study the NeQuick model version 2.0.2 is applied and driven by the daily solar radio flux index F10.7 (see Hochegger et al., 2000;Radicella and Leitinger, 2001;Nava et al., 2008).
The paper is organized as follows: Sect. 2 describes the geostatistical theory of kriging to reconstruct ionospheric parameters, such as STEC and the electron density.In Sect. 3 we outline the chosen validation scenario based on real and simulated ionospheric measurements.For the two chosen periods of quiet and active ionospheric conditions, STEC and differential STEC (dSTEC) are estimated at four independent GNSS ground stations of the International GNSS Service (IGS) and foF2 characteristics are estimated at the two independent ionosonde stations.The results are compared with the background model results and the 2-D TEC maps of the International GNSS Service (IGS) and the DLR's Ionospheric Monitoring and Prediction Center (IMPC).Subse-quently, Sect. 4 comprises our conclusions and an outlook for future activities.

Methodology
Contrary to the commonly applied methods of ionospheric tomography, kriging does not require the voxelization of the ionosphere and hence avoids the calculation of the intersection geometry between voxels and available ray paths.Instead, kriging operates on the space of the measurements and searches for the optimal covariance between the measurements including the estimation of correlation lengths and variance parameters.In the following subsections the ionospheric reconstruction by means of the developed 3-D kriging approaches is explained in more detail.

Kriging of the electron density with a spatial covariance model
In Minkwitz et al. (2015) a detailed description of the tomography approach based on kriging with a spatial covariance model and STEC measurements is given.In the following we briefly summarize this algorithm.Kriging is a best linear unbiased predictor and relies on the estimation of the covariance between the given measurements as well as between the measurements and the parameters of interest.For that purpose a parametric covariance model is necessary.Following the principle behaviour of the ionosphere, a non-stationary and anisotropic covariance model between two electron densities N e (x i ), N e (x j ) at the ECEF coordinates x i , x j ∈ R 3 of the WGS84 reference ellipsoid is chosen as where θ = (θ 1 , θ 2 , θ 3 , θ 4 ) is the unknown covariance parameter vector; E[N e (x i )] and E[N e (x j )] are the expectation values of the electron densities at the coordinates x i and x j ; Cov θ (N e (s i ), N e (s j ))ds j ds i + Cov(ε s i , ε s j ). (3) In order to derive the maximum likelihood estimate (MLE) of the covariance parameter vector θ , we assume that the STEC measurements − −− → STEC = (STEC s 1 , . .., STEC s n ) T form a Gaussian random field with the vector of expectation values E − −− → STEC , calculated by the NeQuick model, and the covariance matrix ( θ ) ij := Cov θ (STEC s i , STEC s j ) with i, j ∈ {1, . .., n}.Neglecting constants, the MLE of θ is given by We remark at this point that the assumption of a Gaussian distribution might lead to negative electron density estimates, especially in higher altitudes with low electron density values.There are different suggestions to suppress these negative estimates, for instance by the inclusion of pseudo-measurements (see Norberg et al., 2016).Alternatively, one might apply probability density distributions that inherently ensure positive electron density estimates, e.g. the log-normal distribution.However, these methods come along with the need for additional linearization steps (see, e.g., Fridman et al., 2006).Based on the optimized covariance parameter θ the electron density N e (x 0 ) at coordinate x 0 is estimated as where E[N e (x 0 )] is the mean of the electron density provided by the NeQuick model, λ is the optimal weight vector and σ x 0 is the covariance vector between the electron density of interest N e (x 0 ) and the measurements.The corresponding estimation error is calculated as Similarly, the STEC s 0 over a given ray path s 0 is given as

Kriging of the electron density with a spatial-temporal covariance model
The extension of the spatial covariance model to a spatialtemporal covariance model aims at the possibility of including not only measurements of the actual reconstruction time t 0 to estimate the electron density N e (x 0 , t 0 ) but also measurements before and after t 0 , e.g.within a sliding window In order to define a spatial-temporal covariance Cov θ (N e (x i , t i ), N e (x j , t j )) it is often assumed that the covariance may be separated into a purely spatial and purely temporal part (see Gneiting et al., 2006, andGneiting andGuttorp, 2010).Based on this assumption and Eq. ( 1) we define the spatial-temporal covariance as with the temporal correlation C t (h t ; θ 5 ), driven by the absolute time difference h t = |t i −t j | and the temporal correlation length θ 5 .Following the suggestions of Kutiev et al. (1999), Muhtarov and Kutiev (1999), Gizawy (2003) and Bust and Mitchell (2008) the temporal correlation is modelled here by an exponential function.Similar to Eqs. ( 4) and ( 5) the unknown parameter vector θ = (θ 1 , θ 2 , θ 3 , θ 4 , θ 5 ) is obtained by MLE.The electron density N e (x 0 , t 0 ) at coordinate x 0 and time t 0 is estimated by the STEC measurements Analogous to Eqs. ( 7) and ( 9) the STEC over a given ray path s 0 at time t 0 can be estimated.We remark at this point that the given framework of a spatial-temporal covariance model provides the opportunity to predict ionospheric key parameters for a time t 0 which is beyond the currently available measurement epochs.However, this capability is not analysed in the scope of this work.

Gradient-enhanced kriging of the electron density
In general, the developed kriging approaches (see Sects.2.1 and 2.2) are able to assimilate direct and indirect linear measurements of the electron density.The peak electron density characteristics derived from ionosonde stations provide both measurement types.The peak density of an ionospheric 1002 D. Minkwitz et al.: Gradient-enhanced kriging of the electron density layer provides us a direct measurement of the electron density, whereas the corresponding peak height represents an indirect measurement, in particular the information that the directional derivative of the electron density D u N e (x i ) at coordinate x i in height direction u is zero.In order to assimilate such information, the covariance between the electron density and the directional derivative as well as between two directional derivatives, i.e. two peak height measurements D u l N e (x i ), D u m N e (x i ) with l, m ∈ {1, . .., p}, is essential.
According to Näther and Simak (2003), Chilès and Delfiner (2012) and Han et al. (2013) these relations are as follows: Consequently, the existence of the partial derivatives of the covariance models in Eqs. ( 1) and ( 8) is required.However, the exponential correlation model, used so far in the definitions of C h , C v , C t , is not differentiable in the origin (i.e. for x i = x j (see Wackernagel, 2003), and hence has to be replaced appropriately.In this study we choose the Matérn correlation model C Matern (h; ν, a), which is defined as with the set of correlation parameters {θ 2 , . .., θ 5 }, the distance between two measurements h and the modified Bessel function K ν √ 2ν h a and the smoothness parameter ν (Rasmussen and Williams, 2006, chapter 4).The existence of the partial derivatives of C Matern (h; ν, a) can be steered by ν and is related to the existence of the mean square partial derivative of the considered random field.We refer the interested reader for more details to Rasmussen and Williams (2006, chapter 4) and Paciorek (2003, chapter 2).For ν = 1/2 the Matérn correlation model corresponds to the exponential correlation model and for ν → ∞ to the Gaussian correlation model (Chilès and Delfiner, 2012, chapter 2.5.1).Since the estimation of ν from the data is notoriously difficult (see Heaton and Peng, 2013), we conveniently choose ν = 9/2 in this work.This choice leads to the explicit representation By means of Eqs. ( 11) and ( 14) the kriging of the electron density N e (x 0 , t 0 ) with STEC measurements and ionosonde measurements can be formulated.For clarity, in the following equation we only include STEC measurements STEC s 1 , . .., STEC s n and F2 layer characteristics, i.e. the peak densities NmF2(x 1 ), . .., NmF2(x p ) and the corresponding gradient information D u 1 N e (x 1 ), . . ., D u p N e (x p ), of the current time t 0 .Hence, we drop the time index t 0 and get The inverse of the covariance matrix between the measurements −1 θ,GEK is composed of block matrices defining the covariances between measurements of the same type as well as of different types.The work flow of GEK with STEC and ionosonde measurements based on a spatial covariance function is illustrated in Fig. 1.In this paper the implementation of the GEK is restricted to the spatial covariance model.

Validation with real and simulated data
For convenience we denote the previously described methods in the following as M I (kriging with spatial covariance), M II (kriging with spatial-temporal covariance) and M III (GEK).In this section two comparisons are conducted.On the one hand, the capability of M I and M II to reproduce the STEC measurements, spatial-temporal STEC variations dSTEC, and the F2 layer critical frequency foF2 is investigated for a period of low and a period of high solar activity.On the other hand, it is shown how the derivative information given by the F2 layer peak measurements of the ionosondes can be used in M III to reduce the height bias that is observed in the electron density profiles estimated by STEC data only.

Validation scenarios
This study is conducted in the European region, covering 40-60 • N and 30 • W-30 • E for two periods: day of year (DOY) 019-022 of the year 2011 and DOY 296-299 of the year 2014.In the following we refer to these periods as period A Ann. Geophys., 34, 999-1010, 2016 www.ann-geophys.net/34/999/2016/and B, respectively.The periods are in the current solar cycle 24 with a F10.7 of around 80 flux units and a maximum KP index of 2.7 for period A, indicating quiet ionospheric conditions.Conversely, period B is characterized by an average F10.7 of around 218 flux units and a maximum KP of 3.7 pointing out an active ionosphere.For both periods the GPS measurements of around 45 IGS stations are acquired and the STEC measurements are calculated according to Jakowski et al. (2011a).These stations are depicted as black and red triangles in Fig. 2. The red-marked IGS stations, i.e. "ieng", "mad2", "pots", and "spt0", serve as validation stations, and their corresponding great circle distance to the closest IGS station, assimilated in the approaches M I and M II, is given in Table 1.The given measurement geometries of the four stations are used as input for M I and M II to estimate the corresponding STEC, but the STEC measurements are not used within the assimilation.In this study, for M II, we choose a 40 min time window with a sampling at the following time steps: {t 0 ± 20 min, t 0 ± 10 min, t 0 ± 5 min, t 0 }.The refinement of the optimal time window is a tricky task and will require further work that is beyond the scope of this paper.
In addition to the validation with STEC, the differential STEC (dSTEC) is calculated as difference of two ionospheric combinations L I (t i ), L I (t E max ): where the phase measurements L 1 and L 2 at the corresponding frequencies f 1 and f 2 are measured along a continuous arc, N L 1 and N L 2 are the ambiguities, λ L 1 and λ L 2 the wavelengths, B SAT and B REC are the satellite and receiver- related differential phase biases, and t E max is the time when the satellite-receiver link has its maximum elevation angle E max during the arc.In this study we consider arcs with E max ≥ 45 • .Since the measurements are taken along a continuous arc, the difference of the ambiguities and differential phase biases can be considered as zero and the small windup effect on the carrier phase can also be neglected.Hence, we get a very precise dSTEC with σ ( ε) less than or of the order of 0.1 TECU.dSTEC contains spatial and temporal changes in the electron density and is a common tool to analyse the performance of different ionospheric reconstructions techniques without influence of the different calibration techniques (see, e.g., Orús-Pérez et al., 2007;Hernández-Pajares et al., 2009;Feltens et al., 2011).The estimates of the absolute STEC and dSTEC obtained by M I and M II at the four validation stations are compared to the corresponding estimates calculated by the background model NeQuick and the IGS and IMPC TEC maps by applying a standard mapping function.For the entire processing an elevation mask of 10 • and a minimum separation of 15 • to E max are applied.
In addition to the validations on STEC level, the critical frequency of the F2 layer foF2 is estimated for both periods by M I and M II working with GPS STEC only.For that purpose the electron density at the ionosonde stations DB049 (Dourbes, Belgium) and RO041 (Rome, Italy) is estimated with an altitude resolution of 10 km (see red hexagons in Fig. 2).The results are compared to the background information and the auto-scaled measurements of DB049 and RO041.For period B the data for DB049 are only partly available, since between DOY 298 and DOY 300 the ionosonde was not fully operational (S.Stankov, personal communication, 2016).
Moreover, for both periods the STEC measurements and F2 layer characteristics of six ionosonde stations are simulated by the Neustrelitz Total Electron Content Model (NTCM; see Jakowski et al., 2011b), the Neustrelitz Peak Density Model (NPDM; see Hoque and Jakowski, 2011) and the Neustrelitz Peak Height Model (NPHM; see Hoque and Jakowski, 2012).These models depend on location, time and the solar radio flux F10.7.Therefore, the models are tuned according to the daily F10.7 and the STEC measurement geometries, and we choose the locations of ionosondes DB049, EB040, JR055, PQ052, RL052 and RL041 as ionosonde coordinates, shown as hexagons in Fig. 2. Based on this simulation environment the principle of M III is verified.For that purpose, the electron density profiles obtained by the background model and M I are compared to M III at the ionosonde station locations of DB049 and RO041 for selected time stamps.

Real data results
Figure 3 illustrates the capability of the different methods to reproduce STEC at the four independent validation stations for period A (top row) and B (bottom row).Besides the comparison between estimated and measured STEC, the scatter plots contain information about the mean error µ, the mean absolute error |µ|, the root mean square error (RMSE) and the 90 % quantile (Q90) in TECU as well as the mean of the absolute relative error |µ % | in percentage.In general, the different solar and geomagnetic conditions of both periods are clearly indicated by the magnitude of the STEC values observed during these periods.Considering the scatter plot of the NeQuick model, it is visible that there is a discrepancy between the measured STEC and estimated STEC for both periods with µ up to −19.4 TECU.In both periods, the bias is reduced by M I and M II resulting into mean values of about 0 TECU.Furthermore, for M II working with STEC measurements of a sliding window, an additional improvement in the error statistics is achieved in relation to M I.For the STEC estimated by the TEC maps of the IMPC, a µ close to zero is also observed.However, for both periods it is observable that M I and M II outperform the IMPC TEC maps.The STEC values calculated by the IGS TEC maps reveal an offset to the measured STEC and hence the error statistics are degraded.At this point, we remark that the phrase "measured STEC" might be misleading since there is no direct measure of the absolute/calibrated STEC along a GNSS ray path.Instead, a calibration procedure is necessary to estimate the satellite and receiver related biases and subsequently the absolute STEC, commonly referred to as STEC measurement.Since the calibration techniques may vary from research group to research group, the level of the TEC maps may also differ and induce offsets between TEC maps of different research groups.This problem is well known and evokes the demand on additional validation methods, for instance based on dSTEC.
Figure 4 displays the RMS of the absolute STEC estimates for all methods at each single IGS validation station.Note the change in the scale of the y axis between the sub-figures in order to keep the legibility.According to Table 1 the validation station "mad2" has the smallest separation to the next IGS station, whose data are used in the assimilation procedure, and "ieng" exhibits the largest distance.It is observable that for M I and M II the RMS at "mad2" is larger than at "ieng".This fact does not contradict the theory, since the   1).The left-hand figure illustrates the period 019-022/2011 and the right-hand figure the period 296-299/2014.The IGS validation stations are (beginning with the smallest great circle distance) "mad2", "spt0", "pots" and "ieng".The STEC estimates are obtained by the NeQuick model (blue), the kriging with spatial covariance M I (green), the kriging with spatial-temporal covariance M II (red), the IMPC TEC maps (cyan) and the IGS TEC maps (magenta).accuracy of the kriging estimator does not only depend on the separation to the next data point.Much more, it depends on the distribution of the ground-station network itself and the consistency of the data.Furthermore, in both periods M II performs best, followed by M I and the RMS for both methods is similar in both periods except for the station "mad2".Table 2 summarizes the mean error µ, the mean absolute error |µ|, the RMS and Q90 of the measured dSTEC minus the estimated dSTEC in TECU for both periods.Considering the |µ|, RMS and Q90 for these periods, all methods perform better than the NeQuick model working solely without assimilation of STEC measurements.According to the validation results of the absolute STEC, the improvement in the estimates by the usage of a spatial-temporal covariance (M II) instead of a spatial covariance (M I) in the kriging is confirmed.For period B, M II shows the best agreement with the measured dSTEC values, followed by M I, IMPC and IGS.Similarly, for period A the kriging with spatial-temporal covariance (M II) performs best, but regarding the RMS and Q90, the IGS results exhibit a lower error than M I.The capability of M I and M II to reproduce spatial-temporal STEC gradients dSTEC will be further analysed in future.In this context the presented approach M II can be improved by taking into account the available gradient information within the covariance as well as in the estimation itself.
Figure 5 shows the estimated foF2 in MHz for both periods obtained by M I (upper row) and M II (lower row) vs. the modelled and the measured foF2 of DB049 and RO041 and the corresponding error statistics.The results reveal the potential improvement in the background model information by the assimilation of STEC measurements with a mean absolute error |µ| and RMSE that is more than halved during period B. When the foF2 estimates of M I and M II are compared, a reduction of the error level is visible at RO041, whereas at DB049 the estimation errors are almost equal for period A and even increased for period B. However, during both periods the estimation error of M I and M II is within the 95 % percent confidence interval of the estimate itself (see grey shaded area (±2σ ) in Fig. 5).The confidence interval is calculated by Eq. ( 6) and can be interpreted as an indicator for the uncertainty in the estimates resulting from the available measurement geometry and the estimated covariance parameters.

Simulation results
In this paper, the results for the F2 layer peak height hmF2 are not shown for M I and M II, since no additional information is included that might improve the bias in the peak height noticed in the analysis of M I in Minkwitz et al. (2015, Fig. 7).In order to tackle this issue, M III is developed and tested for simulation scenarios on DOY 022/2011 at 00:00, 06:00, 12:00 and 18:00 UTC (see Fig. 6).At both ionosonde stations, the simulated F2 layer peak characteristic (black circle) and the simulated STEC measurements are assimilated.The estimated electron density profiles reproduce the characteristics as a local or global extremum.At 06:00 and 12:00 UTC the algorithm M III (red) works quite well, with a correctly shifted hmF2 and NmF2 in comparison to M I (green).However, in particular at 00:00 UTC it is visible that the simulated peak is met but additionally a second higher peak appears in the M III profile.
This behaviour is caused by a combination of the inconsistencies between the STEC and the ionosonde measurements and the already existing bias in the background electron density profile.Inconsistencies between STEC and ionosonde measurements can lead to cases in which the assimilation of ionosonde data alone provides more accurate foF2 results than by the ingestion of STEC and foF2 simultaneously (see McNamara et al., 2011).Figure 7 depicts the estimated electron density profiles, with M III assimilating only ionosonde measurements (red).Taking into account Fig. 6 at 00:00 UTC it becomes clear that the assimilation of the STEC measurements leads to a second peak density above the simulated peak density height.Consequently, the assimilation of both measurements will require non-trivial measurement error models.Furthermore, if the background has a bias, the estimation of the electron density by Eqs. ( 5) and (9) becomes complicated, since it is assumed that the background model represents the expectation value of the electron density.This underlying assumption is used in the majority of the data assimilation approaches relying on a background model.Biases in the background model can result, for instance, from differences in the ionization level between STEC measurements and background model, or offsets between measured and modelled hmF2.Possible solutions for that issue might be preconditioning of the NeQuick model (or general background model) by means of the mea-  sured STEC or the F2 layer characteristics (e.g.Nava et al., 2006;Brunini et al., 2011) or the combination of different, independent background models (e.g.Elvidge et al., 2016).
The definition of an advanced measurement error model and the preconditioning of the background model will be addressed in future work.Moreover, a natural extension of M III will be the ingestion of further peak characteristics that can be assimilated analogously by means of Eqs. ( 15) and ( 16).

Conclusions
In this study the tomography of the ionosphere by means of kriging is presented.Within this scope the kriging with a spatial covariance and STEC measurements (M I), as presented in Minkwitz et al. (2015), is extended to a kriging with spatial-temporal covariance (M II).Both approaches are mainly based on the estimation of optimal spatial and tem-poral correlation lengths between the given measurements and avoid the voxelization of the ionosphere and thus the corresponding calculation of the intersection geometry.Furthermore, the well-known problem of accurate F2 layer peak height estimation by STEC only is addressed.For that purpose, a novel gradient enhanced kriging (GEK, M III) of the ionosphere is developed.The GEK assimilates STEC measurements as well as direct electron density measurements and derivative information, given, for example, by ionosonde measurements (foF2, hmF2).As initial guess, the methods can use an arbitrary 3-D electron density model.In this study, the NeQuick model is applied.
The methods (M I) and (M II) are validated for two periods (DOY 019-022/2011 and 296-299/2014) of different ionospheric conditions in the European region regarding their capabilities to reproduce absolute STEC, spatial-temporal STEC gradients (dSTEC) and the critical frequency of the F2 layer.The validations on STEC level show that (M I) and (M II) crucially improve the initial guess of the background model and that both methods outperform the STEC given by the 2-D TEC maps of the IMPC and IGS.These promising results are confirmed by the dSTEC validation.In addition, the validations regarding STEC and dSTEC indicate that the estimation errors are further reduced by taking into account STEC measurements of a sliding window in (M II) instead of STEC measurements of the current reconstruction time stamp only (M I).Similarly, the residuals between the measured and estimated foF2 characteristics at DB049 and RO041 are clearly reduced by means of M I and M II compared to the initial guess of the NeQuick model.
Furthermore, the estimates of electron density profiles by GEK in a simulation environment show that direct electron density measurements and gradient information can contribute to overcoming deficiencies in the peak height estimation by ground-based STEC measurements only.These first results also underline the necessity of the best possible initial guess, motivating the improvement in the initial guess, for example, by the tuning of the background model with actual F2 layer characteristics.For that purpose, the additional assimilation of radio occultation measurements might help to provide information in regions where the station infrastructure is sparse.In general, their integration is already possible by M III.However, our upcoming work will be focussed on the step from synthetic to real ionosonde data, which will require advanced measurement error models that are able to cope with the inconsistencies between ionosonde measurements themselves and between the ionosondes and the STEC measurements.

Figure 1 .
Figure 1.Work flow of gradient enhanced kriging with a spatial covariance model.

Figure 2 .
Figure 2. Typical measurement scenario over Europe for DOY 296/2014 using the IGS ground-station GNSS network and six ionosondes within Europe: IGS ground stations are depicted as triangles and ionosondes as hexagons.Stations in red are used for validation.Ionospheric piercing points of the GPS-based STEC measurements are shown as blue circles.

Figure 4 .
Figure 4. Root mean square error of the STEC estimation vs. the great circle distance of the IGS validation station to the closest IGS station in the surrounding of the validation station (seeTable 1).The left-hand figure illustrates the period 019-022/2011 and the right-hand figure the period 296-299/2014.The IGS validation stations are (beginning with the smallest great circle distance) "mad2", "spt0", "pots" and "ieng".The STEC estimates are obtained by the NeQuick model (blue), the kriging with spatial covariance M I (green), the kriging with spatial-temporal covariance M II (red), the IMPC TEC maps (cyan) and the IGS TEC maps (magenta).
Figure 5.Comparison of foF2 estimated by the NeQuick model (blue), the kriging with spatial covariance M I (green, upper row) and the kriging with spatial-temporal covariance M II (green, lower row) at the ionosonde station locations of DB049 and RO041 for both periods.

Figure 6 .
Figure 6.Simulation scenario: comparison of electron density profiles estimated by the NeQuick model (blue), the kriging with spatial covariance M I (green) and the GEK with spatial covariance M III (red).The peak characteristics (black) at the ionosonde station locations of DB049 (top row) and RO041 (bottom row) on DOY 022/2011 are simulated by the models NPHM and NPDM.The assimilated STEC is simulated by NTCM.

Figure 7 .
Figure 7. Simulation scenario: comparison of electron density profiles estimated by the NeQuick model (blue), the kriging with spatial covariance M I (green) and the GEK with spatial covariance M III (red) assimilating only with F2 layer measurements at the ionosonde station location of DB049 on DOY 022/2011.
)] represents the (co-)variance; and C h (h h ; θ 2 , θ 3 ) and C v (h v ; θ 4 ) are the horizontal and vertical correlations with the correlation length parameters in the longitude and latitude direction θ 2 , θ 3 and in the vertical direction θ 4 over the WGS84 reference ellipsoid.We use the exponential correlation model for C h (h h ; θ 2 , θ 3 ) and C v (h v ; θ 4 )

Table 1 .
Overview on the validation stations used and their closest IGS station in the surrounding area.The great circle distance between the corresponding two stations is given under the assumption of a spherical Earth with a radius of 6371 km.