© S. Zhang et al., Published by EDP Sciences 2026

1 Introduction

The first generation of the Constellation Observing System for Meteorology, Ionosphere and Climate (COSMIC) had great success as the first multi-satellite radio occultation mission for probing the Earth’s ionosphere with unprecedented global observation density. The COSMIC ionospheric radio occultation (IRO) data have been widely used for ionospheric studies and have become one of the most important sources of observation for data assimilation due to its global coverage (Zhang et al., 2023). Errors in the IRO observations come from assumptions in the inversion algorithm, carrier phase measurement errors, and orbit errors of GNSS and LEO satellites (Schreiner et al., 2007). The quality of COSMIC ionospheric observations was evaluated by different research teams with direct and statistical comparisons against benchmark instruments – ground-based stations of ionosondes and incoherent scatter radars, which are considered as “truth” (Liu et al., 2010; Krankowski et al., 2011; Hu et al., 2014). These comparisons neglect the errors in the benchmark observations.

For the purposes of data assimilation, the specified observation error includes not only the measurement error but also the representativeness error. The representativeness errors are related to the ability of the model to capture the scale-sizes present in the observations, including horizontal and vertical resolutions as well as errors due to the observation forward operator (Janjić et al., 2018). The correct specification of observation errors is crucial in data assimilation systems to avoid over- or under-weighting the observational information. There is a need to compare different observations with models and estimate the errors in the respective data sets. Historically, the assimilation of IRO electron density profile observations has used only theoretical inversion errors and an assumed constant instrumentation error (Yue et al., 2011). Therefore, they do not completely capture the error that needs to be specified for data assimilation, which has two components: representativeness error and intrinsic (or measurement) error.

The three-cornered hat (3CH) method is a statistical approach for estimating random errors for three data sets simultaneously. The 3CH method was originally developed as an N-cornered hat method and was used for the error variance of N atomic clocks (Gray & Allan, 1974). The study assumed that the errors between the tested oscillators were not correlated and that the error variance of one oscillator could be evaluated even if the errors of two other oscillators were large. By the 1990s, the 3CH method was widely used in physics (Premoli & Tavella, 1993; Ekstrom & Koppang, 2006; Griggs et al., 2014).

The 3CH method requires a minimum of three datasets and can be used to estimate the error variance for each of the three datasets. The overriding assumption of this approach is that the errors of the three observing systems are uncorrelated. The 3CH method is an effective tool for accurately estimating the random error variances of three datasets with relatively large sample sizes, similar scales of representation, similar magnitudes of error between them, and low correlations among the data (Rieckh & Anthes, 2018; Sjoberg et al., 2021). With the advent of multiple overlapping datasets and the elevated need for optimal a priori error estimation in recent years, variations and enhancements of the 3CH method have also been used in geoscience observations, data assimilation, and modeled data (O’Carroll et al., 2008; Long et al., 2014; Awange et al., 2016; Anthes & Rieckh, 2018; Liang et al., 2024).

This work provides a better characterization of the COSMIC electron density profile errors than a one-to-one comparison with other observations, with a specific focus on understanding the intrinsic and representativeness errors in order to improve the error specification for data assimilation. Better characterization of the errors should, ideally, improve the impact of assimilating the observations. In this study, five foF₂ datasets, including COSMIC, ionosonde, IRI, NeQuick, and WACCM-X, are used to evaluate their respective error variances. The 3CH method has never been applied to the ionosphere before, and this work provides a new way for estimating ionospheric observation errors. Furthermore, this method can also be applied to the full vertical profile and other observations.

2 Method

In the 3CH method, the dataset X is assumed to be equal to:

$$X=\mathrm{True}+\beta +\epsilon ,$$(1)

where “True” is the true value, ε is the random error, and β is the mean bias relative to the true value.

Suppose that there are three datasets X, Y, and Z. The variances of the differences between the datasets can be written as

$$\mathrm{VAR}\left(X-Y\right)=\mathrm{VAR}\left({\epsilon }_X\right)+\mathrm{VAR}\left({\epsilon }_Y\right)-2\mathrm{COV}\left({\epsilon }_X,{\epsilon }_Y\right),$$(2.1)

$$\mathrm{VAR}\left(X-Z\right)=\mathrm{VAR}\left({\epsilon }_X\right)+\mathrm{VAR}\left({\epsilon }_Z\right)-2\mathrm{COV}\left({\epsilon }_X,{\epsilon }_Z\right),$$(2.2)

$$\mathrm{VAR}\left(Y-Z\right)=\mathrm{VAR}\left({\epsilon }_Y\right)+\mathrm{VAR}\left({\epsilon }_Z\right)-2\mathrm{COV}\left({\epsilon }_Y,{\epsilon }_Z\right).$$(2.3)

Then we can get (see Appendix A in Anthes & Rieckh, 2018):

$$\begin{array}{c}\mathrm{VAR}\left({\epsilon }_X\right)=\frac{1}{2}\mathrm{VAR}\left(X-Y\right)+\frac{1}{2}\mathrm{VAR}\left(X-Z\right)-\frac{1}{2}\mathrm{VAR}\left(Y-Z\right)\\ +\mathrm{COV}\left({\epsilon }_X,{\epsilon }_Y\right)+\mathrm{COV}\left({\epsilon }_X,{\epsilon }_Z\right)-\mathrm{COV}\left({\epsilon }_Y,{\epsilon }_Z\right),\end{array}$$(3.1)

$$\begin{array}{c}\mathrm{VAR}\left({\epsilon }_Y\right)=\frac{1}{2}\mathrm{VAR}\left(X-Y\right)+\frac{1}{2}\mathrm{VAR}\left(Z-Y\right)-\frac{1}{2}\mathrm{VAR}\left(X-Z\right)\\ +\mathrm{COV}\left({\epsilon }_X,{\epsilon }_Y\right)+\mathrm{COV}\left({\epsilon }_Y,{\epsilon }_Z\right)-\mathrm{COV}\left({\epsilon }_X,{\epsilon }_Z\right),\end{array}$$(3.2)

$$\begin{array}{c}\mathrm{VAR}\left({\epsilon }_Z\right)=\frac{1}{2}\mathrm{VAR}\left(X-Z\right)+\frac{1}{2}\mathrm{VAR}\left(Z-Y\right)-\frac{1}{2}\mathrm{VAR}\left(X-Y\right)\\ +\mathrm{COV}\left({\epsilon }_X,{\epsilon }_Z\right)+\mathrm{COV}\left({\epsilon }_Y,{\epsilon }_Z\right)-\mathrm{COV}\left({\epsilon }_X,{\epsilon }_Y\right).\end{array}$$(3.3)

If the covariance (COV) terms are negligible (i.e., the error covariance is small compared to the variance difference between the data sets or the terms approximately cancel), then the error variance for each data set is

$$\begin{array}{c}\mathrm{VAR}\left({\epsilon }_X\right)=\frac{1}{2}\mathrm{VAR}\left(X-Y\right)+\frac{1}{2}\mathrm{VAR}\left(X-Z\right)-\frac{1}{2}\mathrm{VAR}\left(Y-Z\right),\end{array}$$(4.1)

$$\begin{array}{c}\mathrm{VAR}\left({\epsilon }_Y\right)=\frac{1}{2}\mathrm{VAR}\left(X-Y\right)+\frac{1}{2}\mathrm{VAR}\left(Z-Y\right)-\frac{1}{2}\mathrm{VAR}\left(X-Z\right),\end{array}$$(4.2)

$$\begin{array}{c}\mathrm{VAR}\left({\epsilon }_Z\right)=\frac{1}{2}\mathrm{VAR}\left(X-Z\right)+\frac{1}{2}\mathrm{VAR}\left(Z-Y\right)-\frac{1}{2}\mathrm{VAR}\left(X-Y\right).\end{array}$$(4.3)

After the variances are calculated, the standard deviation (STD) can be calculated by $\mathrm{STD}=\sqrt{\mathrm{VAR}}$.

The 3CH method gives estimates of the statistics of random errors (precision or uncertainty), not systematic errors (accuracy). The dominant source of inaccuracy in the 3CH error variance estimates is the correlation of errors among two or more of the datasets. For three (or more) observational datasets that sample the same field with the same footprint (scale representation of the field) and at the same time/location, the assumption of uncorrelated errors is considered to be valid. Error correlation may arise due to inclusion of a model that assimilates the observations being used in the 3CH, representativeness differences, collocation errors, and correlations arising by chance due to a small sample size (Rieckh et al., 2021). A small sample size or a large difference between the errors of the three datasets, even for small error correlations, will also reduce the accuracy of the error estimates (Sjoberg et al., 2021).

Through practical application on various real datasets, a deeper understanding of the 3CH method’s implementation is gained. When using reanalysis products that assimilate observational data, neglecting the covariance terms in certain estimates can affect the results. Substantial evidence illustrates that the estimated error statistics contain valuable information that outweighs the noise introduced by ignoring the covariance term and by the limited data sample (Anthes & Rieckh, 2018). Calculations with simulated data suggest that the effect of error correlations is negligible (produces errors in the estimated error standard deviations less than 10%) for correlation coefficients that are smaller than 0.1, assuming the magnitude of the errors of the correlated datasets is similar (Rieckh & Anthes, 2018).

Representativeness is important to consider when analyzing and interpreting 3CH results. For example, when using combinations of a large number of datasets, representativeness errors may cause different clusters of error estimates that are associated with observations having different resolutions. Interpreting the differences of these clusters depends on whether truth is associated with the high- or low-resolution observations (Rieckh et al., 2021). Additionally, when two datasets with similar vertical resolutions form a triad with a dataset with a different vertical resolution, the 3CH method associates the “truth” with the scale of the two datasets with similar resolutions and penalizes the third dataset (Sjoberg et al., 2021).

The 3CH is also used to estimate error statistics or uncertainties of relevance to data assimilation diagnostics (Semane et al., 2022). When the 3CH method uses the observation, background, and analysis for its three corners, the third corner estimate can recover the negative of the analysis error variance, which can be used as a supplement to traditional residual statistics (Todling et al., 2022).

3 Data

Five foF₂ datasets from 2006 to 2019 are used in this paper. The COSMIC and ionosonde data are collocated using collocation criteria of 600 km and 30 min, and are considered to be at the same location and time. Regardless of the collocation, a representativeness error will be introduced, particularly in regions of significant ionospheric structuring such as high latitudes and in the equatorial post-sunset region, where the ionosphere exhibits substantial horizontal gradients and undergoes dynamic changes on time scales much shorter than 30 min. However, overall, this collocation criterion is reasonable and consistent with prior studies (Cherniak et al., 2021).

The foF₂ values from the IRI, NeQuick, and WACCM-X models are interpolated to the time and location of the occultation tangent point. Before the error calculation using the 3CH method, an error Gaussian is fitted between every two of the three data sets, and outliers are removed if they exceed three times the standard deviation. In this paper, the NmF₂ data are transformed into foF₂ for subsequent calculations:

$$\mathrm{fo}{\mathrm{F}}_2=\sqrt{80.6\mathrm{Nm}{\mathrm{F}}_2}.$$(5)

3.1 COSMIC observations

The FORMOSAT-3/COSMIC mission is a six-satellite IRO constellation dedicated to remotely sensing Earth’s atmosphere and ionosphere, and was launched in 2006 (Schreiner et al., 2007). The occultations yield abundant information the slant total electron content as well as electron density profiles. In this research, we used the inversion result of the COSMIC ionospheric GPS occultation data from the University Corporation for Atmospheric Research (UCAR) COSMIC Data Analysis and Archival Center (CDAAC), which is derived through an Abel inversion. The assumptions and approximations used in the current CDAAC processing include the following: (1) proportionality between refractivity and electron density, (2) straight-line signal propagation, (3) spherical symmetry of electron density, (4) circular satellite orbits, and (5) a first-order estimate of the electron density at the top (Lei et al., 2007).

3.2 Ionosonde observations

Ionosondes are a widely used ionospheric measurement technique. This study uses data from the Digital Ionogram Database (DIDBase) in the Lowell Global Ionosphere Radio Observatory Data Center (LGDC) (Reinisch & Galkin, 2011). The locations of the 100 ionosondes used in this study are shown in Figure 1. The Automatic Real-Time Ionogram Scalar with True height (ARTIST) is an intelligent system developed for the automatic extraction of ionospheric specification data from digisonde ionograms (Galkin & Reinisch, 2008). DIDBase uses ARTIST-5, and ARTIST is the most widely used automatic scaling routine. For 95% of automatically processed ionograms, the value of the foF₂ parameter determined with ARTIST-5 is within the limits of (−0.3,…,+0.3 MHz) from the foF₂ values determined by a human expert (Galkin & Reinisch, 2008). The autoscaling algorithm underestimates the true values of the frequencies at some fixed frequencies. This is caused by gaps in the echo trace that lead to the premature truncation of the automatic layer traces. ARTIST determines foF₂ as the maximum frequency on the scaled O-trace, which is clearly incorrect when that trace is ended prematurely (Stankov et al., 2012). Note that we use auto-scaled ionograms despite their limitations due to requiring a large dataset of observations to use in the 3CH method, which is difficult to obtain using manually scaled ionograms. The use of auto-scaled ionograms and the impact of auto-scaling on the ionosonde errors are important to consider in the present study. The temporal resolution of ionosonde measurements varies among stations, ranging from 1 h to 15 min, and in some cases, 5 min. Ionosonde data are paired with COSMIC observations that occur within a ±30-min window of the ionosonde observations in this paper.

Figure 1 Global distribution of ionosonde sites. Red dots indicate ionosonde site locations. The different colors represent the number of observation data used in this study after pairing with COSMIC IRO.

3.3 IRI model

IRI is a joint undertaking by the Committee on Space Research (COSPAR) and the International Union of Radio Science (URSI) to develop and improve an international standard for the parameters in Earth’s ionosphere (Bilitza et al., 2017). The 2020 version of IRI was used in this study (Bilitza et al., 2022). The core model describes monthly averages of the electron density, electron temperature, ion temperature, and ion composition globally in the altitude range from 60 to 2,000 km. IRI models the ionospheric vertical electron density profile by dividing it into six sub-regions. The boundaries of these subsections are identified by several characteristic points, including the F₂, F1, and E-layer peaks. The foF₂ is based on the long data record accumulated by the worldwide network of ionosondes consisting of 150–170 stations and going back to the fifties. In recent years, the steadily increasing network of digisondes, which provide real-time values of foF₂, has greatly expanded the global coverage for these important ionospheric characteristics. The IRI model is the most popular empirical ionospheric model and is comparable to the best theoretical and assimilation models in predicting foF₂ (Shim et al., 2023).

3.4 NeQuick model

NeQuick is a three-dimensional and time-dependent ionospheric electron density model developed at the Aeronomy and Radio propagation Laboratory of the Abdus Salam International Centre for Theoretical Physics, Trieste, Italy, and at the Institute for Geophysics, Astrophysics and Meteorology of the University of Graz, Austria (Nava et al., 2008). NeQuick-2 was used in this study (Nava et al., 2008). It is a quick-run model, particularly tailored for trans-ionospheric applications, that allows one to calculate the electron concentration at any given location in the ionosphere. Three profile anchor points are used: the E layer peak, the F1 peak, and the F₂ peak, which are modeled in terms of the ionosonde parameters foE, foF1, foF_2, and M(2000)F₂. Each part of the profile is modeled using the Epstein layer formalism (Coïsson et al., 2008). Versions of the IRI model from 2007 and later have included the topside formulation of NeQuick, and have adopted it as the most mature of the different proposals (Bilitza et al., 2022). To compute median values of foF₂ and M(3000)F₂, the NeQuick model uses the ITU-R (formerly called CCIR, 1967) coefficients (Jones & Gallet, 1962).

3.5 WACCM-X model

WACCM-X is a configuration of the NSF NCAR Community Earth System Model (CESM) that extends the atmospheric component into the thermosphere, with a model top boundary between 500 and 700 km, depending on solar activity. It self-consistently resolves the lower and upper atmospheric processes, enabling a more realistic simulation of the climatology and variability (Liu et al., 2018). Physical processes represented in WACCM-X build upon those in regular WACCM, which has a model top at ∼140 km, and in turn is built upon the Community Atmosphere Model (CAM), which goes up to ∼40 km. The ionospheric electrodynamics in WACCM-X are adapted from the TIE-GCM. The ionospheric component in WACCM-X includes modules of the ionospheric wind dynamo, F-region O⁺ transport, and electron and ion temperatures, which are used to calculate heating of the neutral atmosphere through collisions with thermal electrons and ions. WACCM-X can produce better quiet time foF₂ than TIE-GCM does and capture wave-like small increases in foF₂ (Shim et al., 2023). In the present study, we use Specified Dynamics (SD) WACCM-X simulations. The meteorology in the SD-WACCM-X (hereafter just WACCM-X) simulations is constrained to the NASA MERRA-2 reanalysis up to ~50 km. Time-varying Kp and F10.7 are used to specify the high-latitude forcing (electric potential and aurora) and solar irradiance, respectively.

3.6 Dataset difference analysis

To understand the datasets used in this paper, we first compare the ionosonde observations with the other four datasets. Since ionosonde measurements are generally regarded as the most accurate observations, we calculated the mean and standard deviation of the differences between the other four datasets and the ionosonde observations. The mean and standard deviations of the differences are given in Table 1 for all ionosondes globally. The corresponding histograms of the differences are shown in Figure 2. It can be seen that the main differences in all the datasets come from the STD rather than the mean. Especially for IRO, the mean difference is minimal. Since the mean bias is relatively small, the differences for ionosondes and COSMIC IRO, as well as the different models, can be considered to be dominated by a random component. We note that the STD of the IRI-ionosonde and NeQuick-ionosonde differences is partly reflective of these empirical models not adequately capturing the day-to-day variability. However, for the purposes of data assimilation, this can be considered as errors due to the misrepresentation of physical processes in the background model that need to be considered.

Figure 2 Histograms of the differences between ionosonde observations and each dataset. The upper left is COSMIC IRO, the upper right is IRI, the lower left is NeQuick, and the lower right is WACCMX. The same coordinate scale is used.

4 Simulation study

To first verify if the 3CH method is suitable for IRO validation, we conducted a simulation study. The NeUoG model, the electron density model developed at the University of Graz (Leitinger et al., 1997), is used as the “true” background ionosphere for the simulation. The IRO profiles are simulated using an Abel inversion. The IRI and ionosonde data are generated by adding both bias and random errors to the foF2 values from the NeUoG ionosphere. The bias and random errors added to generate the IRI and ionosonde data are provided in Table 2. This results in three datasets with known errors, which we use to verify the accuracy of the errors calculated using the 3CH method.

We simulate three datasets at the locations of a sample of 1,039 IRO events. The occultations are approximately uniformly distributed in latitude and longitude (Fig. 3). In the simulation, UT time of the NeUoG model is set as 12:00, and F10.7 is set as 150 solar flux units. The month in the model is set as January, April, and July, respectively, to represent different seasons. This means that the orbit data of the observed 1,039 IRO events are used in the simulation for each of the three months, resulting in a total of 3,117 occultations. These simulated IROs are inverted using the Abel inversion, and the foF₂ of the inversion results are validated by the true NeUoG model.

Figure 3 The distribution of the simulated IRO observations. The black dots represent the locations of the tangent points for different occultation events.

The IRO column in Table 2 shows the mean and STD of the differences between the “true” ionosphere and the simulated IRO inversions. The pseudo ionosonde column indicates the specified errors for the simulated ionosonde observation data, i.e., the mean systematic error of the simulated ionosonde is −0.001 MHz and the STD is 0.405 MHz. The values were chosen to be equivalent to the error of the automatic scaling of ionograms (Stankov et al., 2012; Stankov et al., 2023). Similar to the pseudo ionosonde column, the pseudo IRI column gives the specified mean and STD of the error of the pseudo IRI dataset, the mean systematic error of simulated IRI is −0.147 MHz, and the STD is 0.815 MHz.

The 3CH method is used to calculate the STDs of the three datasets, and the results are shown in the last row of Table 2. The values in the second (truth) and third (estimated by the 3CH method) rows are approximately equal. The similarity between the specified standard deviation (STD_True) and the 3CH-derived error (STD_3CH) line demonstrates that the 3CH method can effectively estimate the STD of these datasets. It means that the 3CH method could be used to estimate the error of IRO with two anchor datasets, such as ionosondes and IRI.

5 Results and discussion

Before showing and discussing the interpretation of the results derived from different triplet combinations, several aspects of the observations and models should be discussed, as these influence how the results are interpreted. We would expect no error correlation between the COSMIC and ionosonde observations. The exception is if plasma irregularities are present, which can generate errors in both COSMIC and ionosonde observations. The empirical models (IRI and NeQuick) are developed from observations, especially ionosonde observations, and there may therefore exist error correlations between the two models as well as between the models and ionosonde observations. Although COSMIC is also used to develop IRI (Bilitza et al., 2022) IRO data is used for modeling of hmF₂ instead of foF₂. Therefore, the IRO observation is considered to be independent.

It also should be noted that both the empirical and physics-based (WACCM-X) models have relatively coarse resolution and are not expected to resolve the small spatial and temporal scales that are present in the observations. In contrast, the COSMIC and ionosonde observations likely have relatively similar observational footprints, at least in comparison to the models. The horizontal footprint of the IRO observations, however, is greater than that of the ionosondes.

Based on this information and the characteristics of the 3CH method described above, we would expect the COSMIC 3CH error derived from using COSMIC, ionosonde, and a single model to primarily illustrate the intrinsic error of the observations, which includes both measurement and processing errors. On the other hand, the COSMIC 3CH error derived from using a combination of COSMIC observations and two of the models contains error due to the models and observations having different footprints (i.e., resolved spatial and temporal scales), and thus representing different scales of the ionosphere. This error will therefore include the representativeness error that results from the COSMIC observations resolving different spatial and temporal scales than the models.

5.1 Intrinsic error

The results in Figure 4 show the error STD distribution for the COSMIC, ionosonde, and IRI foF₂ in local time (LT) and geomagnetic latitude (MLat), where the absolute value of the MLat is considered. Since the data are matched using COSMIC and ionosonde, and ionosonde is mainly distributed in the mid-latitudes, the maximum values of the sampled numbers are distributed in the mid-latitudes. The x-axis of each subplot is LT at 1-h intervals, and the y-axis is MLat at 10° intervals. Errors are expressed as STD in MHz. The blank grids (white shading) in the error results indicate STDs that are less than zero, which arise from the aforementioned neglect of the error covariance terms. They may occur where the STDs are small, and thus the error due to neglecting the covariance terms can be large enough to result in a negative estimate of the 3CH STD. To avoid abnormal results caused by small sample sizes within a 1-hour interval, Figure 5 shows the COSMIC STD at a 3-hour interval calculated using a combination of COSMIC, ionosonde, and one model (IRI, NeQuick, or WACCM-X).

Figure 4 The STD and sample number in the triplet COSMIC-Ionosonde-IRI. From left to right, the error of COSMIC IRO, the error of ionosonde, the error of IRI, and the distribution of sample size.

Figure 5 The COSMIC STD in the triplet COSMIC-Ionosonde-Model at 3-h interval. The models from left to right are IRI, NeQuick, and WACCM-X.

The results in Figures 4 and 5 reveal that the largest errors occur at low geomagnetic latitudes and during the afternoon. Larger errors are also seen at high latitudes. The magnitude and spatial distribution of the IRO errors are largely independent of the model, providing confidence that the 3CH error estimates are reliable and representative of the intrinsic observational error. The largest errors occurring at low latitudes during daytime are consistent with the fact that the COSMIC IRO errors are primarily dominated by the assumption of spherical symmetry in the Abel inversion (Yue et al., 2010). The high latitude errors may be due to the occurrence of F-region irregularities and GPS cycle slips. The F-region irregularities and GPS cycle slips primarily occur at night (Yue et al., 2016), consistent with the enhanced COSMIC errors in Figure 4. Increased high latitude nighttime errors were also seen in a comparison between in-situ electron density observations and COSMIC IRO observations (Pedatella et al., 2015). The intrinsic errors are larger in the equatorial ionization anomaly (EIA) region than at the equator. Simulations demonstrate that the Abel retrieval method introduces significant errors in electron densities in the low latitude region. It overestimates electron density in the north and south of the EIA crest region, but underestimates the electron density surrounding the EIA crests (Yue et al., 2010). The 3CH method thus estimates COSMIC errors that are largely consistent with previous results. However, it can provide a more quantitative assessment of the errors.

The ionosonde errors in Figure 4 are larger than the COSMIC IRO errors, which is perhaps a surprising result. For most applications, the accuracy of ionosondes is acceptable, with good reliability when auto-scaling the widely used foF₂. The autoscaling error of ionosonde could be ~0.2–0.4 MHz (Stankov et al., 2012; Stankov et al., 2023). This is comparable to the errors seen in Figure 4 at middle latitudes during nighttime. Based on comparison of foF₂ based on auto-scaling and manual processing of ionosonde observations, the majority of the larger errors in the automatically scaled value of foF₂ were found to consistently occur just after dusk (Bamford et al., 2008). The auto-scaled foF₂ data show pronounced underestimation of the soundings in the dusk hours during high and medium solar activity, but overestimation occurs during the winter nights at high and medium solar activity. These two patterns gradually fade away at lower levels of solar activity, when they are replaced by slight underestimation during the daytime (Stankov et al., 2012). Overall, dusk is the time when ionosondes have the largest observational error due to auto-scaling.

Additionally, when interpreting the ionosonde results in Figure 4, it is important to consider the fact that IRO-based electron density profiles should not be interpreted as classical vertical profiles. The geographical projections of the IRO tangent points at the top and at the bottom of a profile may be separated by up to several hundred kilometers (Lei et al., 2007; Wu et al., 2009; Yue et al., 2010). The IRO observations are also based on inversion of the slant-path TEC, which will tend to average over small-scale features in the ionosphere. The geometry of ground-based ionosonde sounding and space-based IRO measurements is thus different, and this comparison suffers from representativeness differences. IRO can represent a larger detection range, while the ionosonde only represents single-point detection, which means the detection of COSMIC IRO is more representative of the spatial scales of the models than the ionosonde observations. As discussed in Section 2, representativeness differences impact the 3CH estimates through the neglected error covariance terms. The estimates of ionosonde uncertainties and the COSMIC uncertainties do not solely reflect the observation errors, but rather the observation errors plus or minus (respectively) errors due to representativeness differences. This is similar to the forecast error situation for RO and radiosonde in the atmosphere (Kuo et al., 2004; Anthes & Rieckh, 2018). Previous studies found that the errors in estimated radiosonde observations are larger than those in RO observations, partly due to the representative errors of the radiosonde, which provide in situ point measurements, while the model data are larger-scale horizontal averages similar to the RO data.

Table 3 compares the COSMIC errors derived from 3CH (Eq. 4) with the results of Cherniak et al. (2021, hereafter referred to as Ch21), who derived COSMIC-2 errors based on comparison with ionosondes under the assumption of no errors in the ionosonde observations (Cherniak et al., 2021). Note that COSMIC-1 and COSMIC-2 electron density profiles are considered to have similar errors. Daytime is defined as 08-20 LT, and nighttime as 20-08 LT. Low latitudes are below 30° MLat, while areas above this are considered mid-latitudes. The first two columns show the mean and STD obtained for the errors between the ionosonde and COSMIC-2 from Ch21. Columns 3 and 4 show the mean and STD for the direct comparison of COSMIC and ionosonde observations used in this study. Columns 5, 6, and 7 show the results of the COSMIC errors when derived based on the 3CH method using a triplet of COSMIC, ionosonde, and one of the models. The smaller error results in Ch21 are attributed to the use of manually scaled ionograms and a collocation criteria of within 5° and 15 min of the closest IRO and only geomagnetically quiet conditions, whereas in this paper, we use auto-scaled ionograms as well as a less strict collocation criteria to guarantee a sufficient number of collocations. The magnitudes and distribution patterns of the 3CH STD errors exhibit similarities with those derived based on direct comparison of COSMIC observations to ionosondes. For example, the errors are larger at low latitudes than at middle and high latitudes, and smaller at night than during the day, which is consistent with the expected COSMIC inversion errors. This gives us confidence that the COSMIC error derived from the 3CH method is generally reasonable. Some of the differences in mean values between low and mid-latitudes may be due to the coverage of COSMIC-1 and COSMIC-2. Therefore, the inclusion of regions with higher latitudes in the COSMIC-1 data in this study may contribute to this difference.

The impact of potential error correlation between the IRI model (the following analysis also applies to the NeQuick model) and ionosondes on the results should also be considered, i.e., the impact of neglecting the COV(ε_sonde, ε_IRI) term. The errors of IRO and ionosonde estimated here can be expressed as (see in Eqs. 2b and 4b, Rieckh & Anthes, 2018)

$$\mathrm{VAR}({\epsilon }_{\mathrm{IRO}}{)}_{\mathrm{est}}=\mathrm{VAR}\left({\epsilon }_{\mathrm{IRO}}\right)+\mathrm{COV}\left({\epsilon }_{\mathrm{sonde}},{\epsilon }_{\mathrm{IRI}}\right),$$(6)

$$\mathrm{VAR}({\epsilon }_{\mathrm{sonde}}{)}_{\mathrm{est}}=\mathrm{VAR}\left({\epsilon }_{\mathrm{sonde}}\right)-\mathrm{COV}\left({\epsilon }_{\mathrm{sonde}},{\epsilon }_{\mathrm{IRI}}\right).$$(7)

This indicates that correlation between the ionosonde and IRI errors leads to a certain degree of overestimation of the IRO error and underestimation of the ionosonde error and IRI error when the 3CH triplet is comprised of IRO, ionosonde, and IRI. However, the results obtained using WACCMX, which is a purely physical model, also show similar errors. This shows the suitability of the analysis and indicates that the neglected error covariance terms between ionosondes and the empirical models do not significantly influence the results.

5.2 Representativeness error

Figure 6 shows the COSMIC 3CH error derived using a combination of COSMIC observations and two of the models. These errors result from a combination of intrinsic and representativeness errors. The effect of the resolution differences on the 3CH estimates can be seen as beneficial, as the error variance for data assimilation should include representativeness errors that result from mismatched observational and model footprints. When specifying error variances for data assimilation, it is common to artificially increase the observational errors in an assimilation system to account for differences in representativeness, hence the need to obtain accurate estimates of the representativeness error. The estimates in these triplets should not be considered to be an estimate of the intrinsic observation error of IRO; however, by accounting for representativeness errors, they are potentially beneficial for the specification of COSMIC error variances in data assimilation systems. Each dataset affects the representativeness of the triplets, and the resolution of the truth should be something like the mean resolution of the triplet, so we cannot know which model dominates the representativeness. However, the resolution of these three model datasets is comparable, so the portion of the 3CH estimated uncertainties due to representativeness is on the correct order of magnitude.

Figure 6 The COSMIC error standard deviation in the triplet COSMIC-Model1-Model2.

Figure 7 illustrates more clearly the difference in errors obtained from 3CH combinations of COSMIC and two models compared to COSMIC, ionosonde, and a single model. Figure 7 shows that the average error based on COSMIC and two models is 0.2–0.4 MHz larger, which indicates that the representativeness error is ~0.2–0.4 MHz. This additional error should be considered in the specification of error variances when assimilating COSMIC IRO observations.

Figure 7 The average error STD of two combinations and their difference. (Left) The average results of Figure 5. (Middle) The average results of Figure 6. (Right) The difference between the two.

5.3 Error normalization

It is also useful to consider the relative errors in the observations, especially given that the electron densities exhibit large day-night and latitudinal differences. We obtain normalized observation errors using the following equation:

$$\mathrm{erro}{\mathrm{r}}_{\mathrm{norm}\left(\mathrm{Mlat},\mathrm{LT}\right)}=\frac{\mathrm{erro}{\mathrm{r}}_{3\mathrm{CH}\left(\mathrm{Mlat},\mathrm{LT}\right)}}{\mathrm{mea}{{\mathrm{n}}_{\mathrm{value}}}_{\mathrm{obs}\left(\mathrm{Mlat},\mathrm{LT}\right)}}\mathrm{*}100\%.$$(9)

Normalization of triplet combinations containing COSMIC IRO by the above equation yields the results shown in Figure 8. The structure of the relative errors in Figure 8 is clearly different than the absolute errors shown in Figure 4. For the results of the COSMIC-Ionosonde-Model combination, the relative error is small in the low-latitude region during the daytime, while it is larger for the low-latitude region and the high-latitude region during the nighttime. The largest relative error during nighttime is also observed in the COSMIC-Model1-Model2 combination.

Figure 8 COSMIC relative error after normalization. The first row is COSMIC-Sonde-Model, the second row is COSMIC-Model1-Model2.

5.4 Variability in COSMIC IRO errors

5.4.1 Variation of COSMIC errors with season

The main source of error in the COSMIC inversion is the assumption of spherical symmetry in the ionosphere. The effect of this assumption is more drastic when the horizontal gradient of ionospheric electron density variation is large. Therefore, the COSMIC errors are expected to vary seasonally. The following results use June–August as summer and December–February as winter in the Northern Hemisphere, and June–August as winter and December–February as summer in the Southern Hemisphere. March–May and September–November are considered equinox months. The statistical results are obtained based on this seasonal categorization. In Figure 9, the first row of results all show the triplet COSMIC-Ionosonde-Model results, characterizing the COSMIC intrinsic errors primarily, and the second row of results shows the triplet COSMIC-Model1-Model2 results, demonstrating the impact of the COSMIC representativeness error.

Figure 9 COSMIC 3CH error in different seasons. The first group of 6 subplots are for equinox, the second group is for winter, and the third group is for summer.

The intrinsic errors of COSMIC are smallest in the summer, largest in the winter, and in-between during equinox, which is consistent with prior simulation and observation results (Wu et al., 2009). However, the representativeness error is smallest in winter and largest in spring and autumn. This may be related to seasonally dependent errors in the models, which could lead to the representativeness error varying with season. It could also arise due to a seasonal dependency in the ability of the models to capture the relevant spatial scales that are present in the observations. The seasonal dependency of the representativeness errors is relatively large, indicating that the seasonal dependency may need to be considered to optimally assimilate the COSMIC IRO observations.

5.4.2 Variation of COSMIC errors with solar activity

Figure 10 shows the COSMIC errors during high and low solar activity. The results are based on considering 2006–2010 and 2017–2019 as low solar activity years and 2011–2016 as high solar activity years. The overall latitude-local time structure of the results is similar to the previously presented results, demonstrating that the overall structure of the errors is largely independent of solar activity.

Figure 10 COSMIC 3CH error in high and low solar activity years. The first 6 subplots are for high solar activity years, and the bottom 6 subplots are for low solar activity years.

Although the distribution of the errors is similar to the previous results, the errors are larger in the high solar activity years than in the low solar activity years. This is most notable for the intrinsic errors (top rows), which are significantly larger during high solar activity years. This is primarily due to the increase in background densities, and the relative errors do not change significantly with solar activity. The average normalized intrinsic error is about 6.37% in high solar activity years and about 6.44% in low solar activity years. The percentage represents the foF₂ errors relative to the original foF₂ values, which removes the influence of absolute differences caused by varying solar activity. The average normalized results for COSMIC and the combination of the two models are about 12.82% in high solar activity years and about 11.82% in low solar activity years. This shows that the relative error is similar under different solar activity conditions. The spherical symmetry assumption is not fully fulfilled in regions with significant horizontal ionospheric gradients, and the gradients may be stronger during high solar activity, leading to larger errors. Simulations of occultation observations show that both the relative and absolute errors in NmF₂ increase in high solar activity years (Wu et al., 2009). The inconsistency in relative errors may be due to the way of the comparison. The errors in this previous simulation study are from directly comparing IRO inversion results with the background model. The model’s ability to accurately capture the effect of solar activity may have caused this difference. As a consequence, the performance of the Abel inversion is degraded at low altitudes and low latitudes (Pedatella et al., 2015), and it may cause some large-scale pseudo-features, such as the two plasma caves underneath the EIA crest (Yue et al., 2012). Thus, the enhancement of solar activity can lead to increased errors.

6 Conclusion

In this study, we employ the 3CH method to analyze error statistics present in foF₂ datasets of COSMIC IRO, ionosondes, IRI, NeQuick, and WACCM-X. We compare the foF₂ datasets from 2006 to 2019 with collocated COSMIC and ionosondes data at 600 km and 30-min intervals and categorize these events based on LT and MLat.

This paper applies the 3CH method to ionospheric observations and models, exploring its applicability to ionospheric data. We analyzed the limitations in the applicability of the method due to dataset error correlations and differences in representativeness. The consistency of results across different triplet combinations and their alignment with previous COSMIC IRO validation studies support the validity of the results obtained through this method (Liu et al., 2010; Krankowski et al., 2011; Hu et al., 2014). The error analysis of the present study can help to understand the errors that need to be specified to most effectively assimilate the observations.

Based on the different triplet combinations results from the 3CH method, we obtain the following conclusions for COSMIC IRO observations:

1. The error varies with LT and MLat, with all foF₂ error variances being largest in the afternoon in low-latitude regions. This is consistent with the error characteristics of the Abel inversion used to retrieve the COSMIC electron density profiles.

2. The representativeness error of IRO with several common models is structurally similar to the intrinsic error of IRO, ranging from ~0.2 to 0.4 MHz.

3. The intrinsic errors of COSMIC all show a minimum in summer and a maximum in winter, while the representativeness errors are the smallest in winter. The absolute COSMIC errors are larger in high solar activity years than in low solar activity years. This is related to the gradient of ionospheric electron density variation.

The results of the present study lead to a better understanding and application of IRO and demonstrate the applicability of the 3CH method to estimate errors of ionospheric observations. This approach can be extended to full vertical profiles or other observations to provide a new perspective for estimating errors that are particularly useful for deriving errors for use in data assimilation systems.

Acknowledgments

The authors thank the COSMIC and GIRO teams for providing the observational data used in this study. We also thank the developers and maintainers of the IRI, NeQuick, and WACCMX models. The editor thanks two anonymous reviewers for their assistance in evaluating this paper.

Funding

This work was originally supported by Chinese National Science Foundation Grants 12241101, 42174192, and 41204137, and by the Strategic Priority Research Program through Chinese Academy of Sciences Grant XDA17010301. The revision of this manuscript was completed during a postdoctoral appointment of SZ at University of Colorado Boulder.

Data availability statement

The ionosonde data are available from GIRO (https://giro.uml.edu/didbase/scaled.php). The data in this work can be downloaded online at https://doi.org/10.5281/zenodo.12688849. The COSMIC IRO data is available from the University Corporation for Atmospheric Research (UCAR) COSMIC Data Analysis and Archival Center (CDAAC) (https://data.cosmic.ucar.edu/gnss-ro/cosmic1/).

References

All Tables

All Figures

	Figure 1 Global distribution of ionosonde sites. Red dots indicate ionosonde site locations. The different colors represent the number of observation data used in this study after pairing with COSMIC IRO.
In the text

	Figure 2 Histograms of the differences between ionosonde observations and each dataset. The upper left is COSMIC IRO, the upper right is IRI, the lower left is NeQuick, and the lower right is WACCMX. The same coordinate scale is used.
In the text

	Figure 3 The distribution of the simulated IRO observations. The black dots represent the locations of the tangent points for different occultation events.
In the text

	Figure 4 The STD and sample number in the triplet COSMIC-Ionosonde-IRI. From left to right, the error of COSMIC IRO, the error of ionosonde, the error of IRI, and the distribution of sample size.
In the text

	Figure 5 The COSMIC STD in the triplet COSMIC-Ionosonde-Model at 3-h interval. The models from left to right are IRI, NeQuick, and WACCM-X.
In the text

	Figure 6 The COSMIC error standard deviation in the triplet COSMIC-Model1-Model2.
In the text

	Figure 7 The average error STD of two combinations and their difference. (Left) The average results of Figure 5. (Middle) The average results of Figure 6. (Right) The difference between the two.
In the text

	Figure 8 COSMIC relative error after normalization. The first row is COSMIC-Sonde-Model, the second row is COSMIC-Model1-Model2.
In the text

	Figure 9 COSMIC 3CH error in different seasons. The first group of 6 subplots are for equinox, the second group is for winter, and the third group is for summer.
In the text

	Figure 10 COSMIC 3CH error in high and low solar activity years. The first 6 subplots are for high solar activity years, and the bottom 6 subplots are for low solar activity years.
In the text