Methods for validation of random uncertainty estimates and their applications to ozone profiles from limb-viewing satellite instruments

For satellite measurements of atmospheric composition, the random uncertainty estimates provided by retrieval algorithms might be imperfect due to various approximations used in the retrievals or the presence of unknown error sources. This paper presents an overview of the methods used for the validation of random uncertainty estimates. All methods discussed in this study are categorized, and assumptions and limitations of each method are discussed. This overview evaluates these methods in application to ozone profile measurements from limb and occultation satellite instruments and provides practical illustrations of random uncertainty validation.

1 Introduction

In nearly all data analyses, such as data comparisons, aggregating/combining/merging data, data assimilation etc., information about data uncertainty is needed. Such characterization of uncertainty would ideally include both systematic and random components, as well as spatio-temporal resolution of the data, as discussed in von Clarmann et al. (2020). Validation of uncertainty estimates is needed, especially if the measurement uncertainty cannot be fully characterized or is based on assumptions. This is typical for remote-sensing measurements, which use retrievals of atmospheric parameters that solve inverse problems. The random uncertainty (i.e. the component of uncertainty that varies randomly and independently between repeated measurements) of the remote sensing measurements is usually estimated via propagation of instrumental noise and other random uncertainties through the inversion algorithm. These estimates, which are sometimes referred to as “ex-ante” errors (von Clarmann, 2006) (other terms are “prognostic”, “predicted”, “inductive”, or “bottom-up”), can be imperfect due to various approximations used in retrievals, or due to the presence of random components in parameter uncertainties.

The aim of this paper is to provide an overview of the methods for validation of random error component. The terms “error” and “uncertainty” are used interchangeably in our paper. We extend the overview of such methods presented in the introduction of Sofieva et al. (2014) and illustrate them using ozone profile retrievals from Earth-orbiting limb and occultation instruments. In our paper, we discuss the applicability and limitations of each method, with the focus on ozone profile retrievals from satellite measurements.

This paper contributes to the APARC (Atmospheric Processes and their Role in Climate) activity TUNER (Towards UNified Error Reporting) https://aparc-climate.org/activities/tuner/ (last access: 13 February 2026).

2 Data

To illustrate the methods employed for validation of random uncertainties, we use ozone profiles retrieved from several limb and occultation measurements. The summary of the datasets is collected in Table 1, and the principles of uncertainty estimates are described below.

Table 1 Information about the ozone profile datasets used in the paper.

2.1 Atmospheric Chemistry Experiment -Fourier Transform Spectrometer (ACE-FTS)

ACE-FTS is a solar occultation instrument operating on board SCISAT satellite from February 2004 to present. ACE-FTS is a high-spectral resolution Fourier-transform spectrometer in the infrared (2.2–13 µm) wavelength regions. Ozone profile retrievals are performed via the global non-linear fitting of the observed spectra. In retrievals from ACE-FTS, estimated random uncertainties of ozone profiles are the fitting errors from the least-squares inversion process (Boone et al., 2005; Sheese et al., 2022). The mean relative random uncertainties are estimated to be lower than 3 % between 12 and 62 km and typically less than 2 % around 30–35 km. Relative uncertainties are slightly higher in polar regions. The ACE-FTS ozone uncertainty estimates slightly grow with time, going from ∼ 1.7 % in the middle stratosphere in the beginning of the mission to ∼ 2.0 % in the recent period. The vertical resolution of ACE-FTS ozone profiles is estimated to be ∼ 3 km.

2.2 Global Ozone Monitoring by Occultation of Stars (GOMOS)

GOMOS was a stellar occultation instrument operated on board Envisat in 2002–2012 (Bertaux et al., 2010). Ozone profiles are retrieved from UV-VIS spectrometer measurements using two-step inversion (Kyrölä et al., 2010): the spectral inversion followed by the vertical inversion. In the vertical inversion, a Tikhonov-type target-resolution regularization is applied (Kyrölä et al., 2010; Sofieva et al., 2004), so that the vertical resolution of retrieved profiles is the same for all occultations. In this paper we use the GOMOS ozone profiles processed with ALGOM2s v.1 Scientific Processor (Sofieva et al., 2017). The error propagation scheme is similar to that used in GOMOS IPF v.6 processor (Kyrölä et al., 2010; Tamminen et al., 2010), as the ALGOM2s ozone profiles are identical to those of IPF v.6 in the stratosphere and differ only in the UTLS. The error estimates (square roots of the diagonal elements of the covariance matrix) are provided in the Level 2 data. The covariance matrix of retrieved profile uncertainties is obtained via Gaussian error propagation through the GOMOS inversion, see Tamminen et al. (2010) for details. Both noise and the dominating random modelling error (due to scintillations) are taken into account in GOMOS inversion. Thus, error estimates provided in Level 2 files represent the total random uncertainty estimates (if neglecting the random part of the parameter errors, which are expected to be minor compared to the abovementioned random error sources).

The random uncertainties of GOMOS ozone profiles depend on stellar brightness, spectral class and obliquity of occultation. They are typically in the range from 0.5 % to 5 % in the stratosphere. Examples of typical uncertainties of GOMOS ozone profiles can be found in Tamminen et al. (2010). An extensive validation of GOMOS random uncertainty estimates is performed and reported in Sofieva et al. (2014). It was shown that GOMOS random uncertainty estimates are realistic for not-dim stars. Due to instrument ageing, GOMOS random uncertainty estimate grow with time, especially for dim stars (Tamminen et al., 2010).

2.3 Michelson Interferometer for Passive Atmospheric Sounding (MIPAS)

MIPAS was an infrared limb emission spectrometer operated on board Envisat in 2002–2012. In this work, we use MIPAS v8 data retrieved with the IMK/IAA processor. Ozone profile retrievals is performed via constrained inverse modelling of limb radiances (von Clarmann et al., 2009). The random errors provided with version 8 IMK MIPAS data consist of the propagated covariance of the spectra, and further uncertainties of parameters used in the retrievals (like uncertainty of the temperature or the line-of-sight pointing) which are of random nature. A detailed description of how these random errors were calculated is provided by von Clarmann et al. (2022) and Kiefer et al. (2023). According to Kiefer et al. (2023, Supplement), the total random error is up to about a factor of 2 to 3 (even >3 for unfavorable conditions like polar winter) higher than the pure propagated measurement noise error in the lower part of the stratosphere (up to about 30 km), and by a factor of 1.1 to 1.4 higher in the upper part of the stratosphere. For illustrations in our paper, if not specified explicitly, we use the total random uncertainty.

2.4 Microwave Limb Sounder (MLS)

MLS is a microwave limb emission spectrometer operating on board the Aura satellite since 2004. This paper uses the Aura MLS “Version 5” dataset (Livesey et al., 2022), retrieved using the same tomographic retrieval algorithm employed for all previous MLS data versions (Livesey et al., 2006; Read et al., 2006; Schwartz et al., 2006). The Level 2 data products consist of vertical profiles spaced 1.5° along the orbital track, with pressure as the vertical coordinate. Each profile is accompanied by a separate profile reporting the estimated precision uncertainty (i.e., random uncertainty) in the profile. This is based on the square root of the diagonal of the solution covariance matrix from the optimal estimation-based retrieval and thus overestimates the scatter in the geophysical products in cases where the a priori information and other regularization constraints contribute significantly to the results (typically at the upper end of the useful vertical range of each product). Systematic errors due to uncertainties in instrument calibration, spectroscopy, and other parameters have been quantified through multiple perturbation studies detailed in Read et al. (2006), and are documented for each version, along with averaging kernels and rules for data use and screening, in a data quality document (Livesey et al., 2022 for version 5).

2.5 Ozone Mapper and Profile Suite – Limb Profiler (OMPS-LP)

OMPS-LP has been operating on board the Suomi- National Polar-orbiting Partnership satellite since 2012. It measures scattered solar light in the limb-viewing geometry. Ozone profiles are retrieved using measurements at UV and Visible wavelengths. For illustrations in this paper, we use three ozone profile datasets from the OMPS-LP instrument: OMPS-LP USask 2D v1.3.0 processed at the University of Saskatchewan, OMPS-LP UBr v4.1 processed at the University Bremen, and the NASA processor v2.6.

The OMPS-LP USask two-dimensional retrieval process uses Gaussian error propagation to estimate the covariance of the retrieved solution due to measurement noise (ignoring the smoothing error). The reported precision is the square root of the diagonal elements of the converged solution covariance matrix. The measurement noise is assumed to be a constant 1 % at all altitudes and wavelengths.

OMPS-LP UBr v.4.1 Level 2 data provide covariance matrices for the retrieved profiles which are obtained by the propagation of the measurement noise errors. The latter are estimated from spectral fit residuals obtained during the pre-processing step. Contributions of the parameter errors are estimated using the Monte-Carlo approach for a set of representative observations and reported by Arosio et al. (2022) along with total error budget estimations.

The NASA OMPS-LP v2.6 algorithm (Kramarova et al., 2024; Rault and Loughman, 2013) retrieves ozone profiles by employing the second order Tikhonov regularization method. The estimated precision for each profile retrieval is calculated using the square roots of diagonal elements of the solution covariance matrix. Systematic errors related to uncertainties in altitude registration, some algorithmic parameters, and a priori profiles have been evaluated and reported by Kramarova et al. (2024); Moy et al. (2017).

2.6 Optical and Spectroscopic Remote Imaging System (OSIRIS)

OSIRIS has been operating on board the Odin satellite since 2001. It uses measurements of scattered solar light for retrievals of ozone and other constituents' profiles. OSIRIS V7.2 data uses standard Gaussian error propagation to estimate the uncertainties in the retrieved ozone profiles. The covariance matrix is calculated through propagation of the measurement noise that is estimated through counting statistics. The reported precision is the square root of the diagonal elements of the covariance matrix.

2.7 Stratospheric Aerosol and Gas Experiment (SAGE)

SAGE II and SAGE III/ISS are solar occultation instruments. The retrieval algorithms for both SAGE II v7.0 (Damadeo et al., 2013) and SAGE III/ISS v5.3 (see SAGE III ATBD, 2002 and Wang et al., 2020) are similar. Ozone profile retrievals are based on the global fit using visible wavelengths (Chappuis band). The uncertainties are computed from the statistical distribution of observations in the L1 transmission data and then propagated through the retrieval algorithm via Gaussian error propagation into the uncertainties for the L2 profiles of ozone and other species. Typical random retrieval uncertainties are within 1 % between 18 and 52 km.

2.8 Scanning Imaging Absorption Spectrometer for Atmospheric Cartography (SCIAMACHY)

SCIAMACHY was a passive remote sensing spectrometer operated on Envisat over 2002–2012. It measured scattered solar light in the ultraviolet, visible and near infrared wavelength region (240–2380 nm). Ozone profile retrievals are performed via constrained inverse modelling of limb radiances using UV and visible spectra. SCIAMACHY V3.5 Level 2 data, which are used in this work, provide covariance matrices for the retrieved profiles which are obtained by the propagation of the measurement noise errors. The latter are estimated from spectral fit residuals obtained during the pre-processing step. Contributions of the parameter errors are estimated using the Monte-Carlo approach for a set of representative observations. Rahpoe et al. (2013) described the method to estimate the parameter errors and reported resulting values along with the total error budget estimation for the precursor retrieval version (V2.5). The results for V3.5 are expected to be similar.

3 Methods for validation of random uncertainty estimates

3.1 Specifics of satellite measurements of atmospheric composition

In remote sensing, retrieved parameters result from solving the inverse problem. The reported random uncertainties are usually estimated via propagation of instrumental noise and other random errors through the inversion algorithm. Error estimates might be wrong if the source uncertainties used for propagation are not known well enough. If some of the error sources are not characterized and the corresponding uncertainties are not considered, the reported uncertainty is underestimated.

The normalized χ² statistics, ${\mathit{\chi }}_{\mathrm{norm}}^{\mathrm{2}}$, is commonly used for assessing the adequacy of the theoretical description of measurements (forward model) and as an indication of the correctness of random uncertainty estimates. ${\mathit{\chi }}_{\mathrm{norm}}^{\mathrm{2}}$ is usually evaluated as:

$$\begin{array}{}\text{(1)}& {\mathit{\chi }}_{\mathrm{norm}}^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{N-p}}{\left(\mathit{y}-{\mathit{y}}_{\mathrm{mod}}\right)}^T{\mathbf{S}}_{y,\mathrm{random}}^{-\mathrm{1}}\left(\mathit{y}-{\mathit{y}}_{\mathrm{mod}}\right),\end{array}$$

where y is the vector of observed parameters (e.g., spectrally resolved radiance values, transmittances), y_mod is the vector of modelled (theoretical) measurements, S_y,random is the covariance matrix of random measurement errors, N is the number of measurements and p is the number of retrieved parameters (e.g., Bevington and Robinson, 2003; Taylor, 1997). ${\mathit{\chi }}_{\mathrm{norm}}^{\mathrm{2}}$ is also called “ χ² per degree of freedom”. If the theoretical model describes the experimental data correctly and the measurement errors are properly defined, ${\mathit{\chi }}_{\mathrm{norm}}^{\mathrm{2}}\approx \mathrm{1}$, ideally. Very large ${\mathit{\chi }}_{\mathrm{norm}}^{\mathrm{2}}$ values indicate underestimated random uncertainties, while ${\mathit{\chi }}_{\mathrm{norm}}^{\mathrm{2}}$ smaller than 1 imply that they are overestimated (Bevington and Robinson, 2003). This simple analysis of ${\mathit{\chi }}_{\mathrm{norm}}^{\mathrm{2}}$ has helped to discover a missing random error component in early GOMOS retrievals, which was due to uncorrected residual scintillation, and to parameterize it in later processing versions (Sofieva et al., 2010; Tamminen et al., 2010). In the first step of the GOMOS retrievals – the spectral inversion (Kyrölä et al., 2010) – ${\mathit{\chi }}_{\mathrm{norm}}^{\mathrm{2}}$ is evaluated at each tangent altitude independently; y and y_mod in the GOMOS case are measured and modelled transmittance spectra, N is the number of spectral pixels (maximum 1416) and p is the number of fitted parameters (slant column densities for ozone, NO₂, NO₃, and aerosol parameters). Figure 1 compares ${\mathit{\chi }}_{\mathrm{norm}}^{\mathrm{2}}$in the GOMOS retrievals in the set of oblique occultations from the brightest star Sirius when residual scintillation errors are ignored (blue) or considered (red).

Figure 1 Adapted from (Sofieva et al., 2010): ${\mathit{\chi }}_{\mathrm{norm}}^{\mathrm{2}}$ in GOMOS retrievals from a set of oblique occultations of Sirius (August 2003, 66° S, obliquity angle is ∼ 25°). Blue: random error due to residual scintillation is ignored, red: modelling errors are taken into account in the retrievals. Dots: values in individual occultations, bold lines indicate median values for the sets. The black dashed line indicates ${\mathit{\chi }}_{\mathrm{norm}}^{\mathrm{2}}=\mathrm{1}$.

A similar analyses of ${\mathit{\chi }}_{\mathrm{norm}}^{\mathrm{2}}$ applied to GOMOS IPF v6 data has spotted ${\mathit{\chi }}_{\mathrm{norm}}^{\mathrm{2}}<\mathrm{1}$ at upper altitudes in case of dim stars, and a more detailed analysis identified the reason and the influence on uncertainty estimates (Sofieva et al., 2014).

In some retrievals, the uncertainty estimate is derived from the fit residuals, as it is done for SCIAMACHY and OMPS-LP UBr ozone retrievals (e.g., Arosio et al., 2022). Such an approach forces ${\mathit{\chi }}_{\mathrm{norm}}^{\mathrm{2}}$ to be close to 1 and can provide an estimate of the uncertainties of measurements in the case one does not trust the information on measurement uncertainty contained in the Level-1 data.

${\mathit{\chi }}_{\mathrm{norm}}^{\mathrm{2}}$ is a statistical and integral characteristic that indicates not only correctness of the random uncertainty estimates but also consistency of measurements with the forward model used for inversion of atmospheric parameters, which is itself often incomplete, or reflects inaccurate knowledge of instrument calibration and/or spectroscopic parameters. In some cases, a priori/regularization terms are included into ${\mathit{\chi }}_{\mathrm{norm}}^{\mathrm{2}}$ too. Sometimes, this metric is used to determine how strong the regularization should be. Therefore, analyses of ${\mathit{\chi }}_{\mathrm{norm}}^{\mathrm{2}}$ at the measurement level are useful, but they cannot substitute validation of uncertainty estimates of the retrieved ozone profiles.

The retrieved ozone profiles are characterized not only by uncertainty estimates but also by vertical resolution and, in the case of tomographic retrievals, by horizontal resolution. This means that, for such retrievals, the impacts of radiance noise on adjacent data points are not independent, and the proper characterization of associated uncertainties is obtained using their covariance matrix.

If the retrieval is performed with the Bayesian maximum a posteriori estimates (Rodgers, 2000; von Clarmann, 2006), a data correlation can also arise due to the usage of a priori information. These aspects should be taken into account when validating uncertainties.

3.2 General strategy

In the laboratory, experimental precision estimates can be obtained using repeated measurements under the same conditions: the sample variance $s^{\mathrm{2}}=\mathrm{var}\left(x\right)=〈(x-〈x〉{)}^{\mathrm{2}}〉$ (angular brackets denote the mean hereafter) approaches the variance of random error distribution σ² (i.e., squared precision) when the size of sample N tends to infinity. For different samples of size N, the values of sample variance will vary due to different random error realization. An ensemble of sample variances is a random variable with a distribution depending on noise variance σ² and N. The quantity $\frac{(N-\mathrm{1})s^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}$ has a χ² distribution with N−1 degrees of freedom, ${\mathit{\chi }}_{N-\mathrm{1}}^{\mathrm{2}}$ (e.g., Bevington and Robinson, 2003; Taylor, 1997). For large N, ${\mathit{\chi }}_{N-\mathrm{1}}^{\mathrm{2}}$ distribution can be approximated by a Gaussian distribution with variance 2N thus var$\left(\frac{(N-\mathrm{1})s^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}\right)=\mathrm{2}N$. This gives the uncertainty of the experimentally estimated random error

$$\begin{array}{}\text{(2)}& \mathrm{var}\left(s^{\mathrm{2}}\right)\approx {\mathit{\sigma }}^{\mathrm{4}}{\displaystyle \frac{\mathrm{2}}{N}}.\end{array}$$

In contrast with many laboratory experiments, geophysical observation conditions cannot be kept exactly constant for atmospheric measurements. Therefore, the sample variance contains a contribution from the natural variability ${\mathit{\sigma }}_{\mathrm{nat}}^{\mathrm{2}}$:

$$\begin{array}{}\text{(3)}& s^{\mathrm{2}}\approx {\mathit{\sigma }}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{nat}}^{\mathrm{2}}.\end{array}$$

For validation of uncertainty estimates, ${\mathit{\sigma }}_{\mathrm{nat}}^{\mathrm{2}}$ should be minimized by selecting collocated measurements or it should be estimated from independent sources (for example, from a chemistry-transport model, CTM).

Approaches for validation of error estimates usually rely on the variance of the difference, $s_{\mathrm{12}}^{\mathrm{2}}=\mathrm{var}(x_{\mathrm{1}}-x_{\mathrm{2}})$, in a set of collocated measurements x₁ and x₂:

$$\begin{array}{}\text{(4)}& s_{\mathrm{12}}^{\mathrm{2}}={\mathit{\sigma }}_{\mathrm{0},\mathrm{nat}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}.\end{array}$$

In Eq. (4), ${\mathit{\sigma }}_{\mathrm{0},\mathrm{nat}}^{\mathrm{2}}$ stands for the natural variability within a space-time collocation window (note that ${\mathit{\sigma }}_{\mathrm{0},\mathrm{nat}}^{\mathrm{2}}$, which represents the mismatch uncertainty, is different from ${\mathit{\sigma }}_{\mathrm{nat}}^{\mathrm{2}}$ in Eq. (3) that represent natural variability in a certain location).

It is important to note that for the vertically resolved ozone profile data involved, calculating differences and combining data require harmonization of data representations in terms of physical quantities and vertical sampling at least. As the satellite data result from a retrieval process, knowledge of prior information and averaging kernel matrices in principle allows retrieval differences to be accounted for as well. Keppens et al. (2019) provide an overview of harmonization operations for atmospheric profile observations, covering vertical representation matching, vertical smoothing matching, and retrieval matching (essentially the prior information contributions). The effect of these manipulations on the information content and uncertainty budget of the original data is extensively discussed in that work and will not be repeated here. In the following, we assume that all profiles are presented in a similar vertical resolution, so vertical smoothing difference errors between profiles can be neglected. For the illustrations in this paper, the harmonization of the vertical resolution is not needed, as the satellite limb profiles considered in this study have similar vertical resolution, see Table 1.

We divide the methods for random uncertainty validation into two groups depending on what kind of data are used: (1) from the same instrument and (2) from different instruments.

3.3 Using collocated measurements from the same instrument

For perfectly collocated measurements (${\mathit{\sigma }}_{\mathrm{0},\mathrm{nat}}^{\mathrm{2}}\approx \mathrm{0}$) from the same instrument with the same precisions ${\mathit{\sigma }}_{\mathrm{1}}={\mathit{\sigma }}_{\mathrm{2}}=\mathit{\sigma }$, Eq. (4) is reduced to $s_{\mathrm{12}}^{\mathrm{2}}\approx \mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}$, thus allowing validation of the uncertainty estimate ${\widehat{\mathit{\sigma }}}^{\mathrm{2}}=s_{\mathrm{12}}^{\mathrm{2}}/\mathrm{2}$. In this estimate, random errors in x₁ and x₂ are assumed to be uncorrelated. This uncertainty validation method was realized, for example, for closely collocated MIPAS ozone profiles (Piccolo and Dudhia, 2007) and OSIRIS ozone measurements (Bourassa et al., 2012). The uncertainty of this experimental precision estimate is defined by the uncertainty of sample variance $s_{\mathrm{12}}^{\mathrm{2}}$.

There are several limitations associated with this method. First, natural variability ${\mathit{\sigma }}_{\mathrm{0},\mathrm{nat}}^{\mathrm{2}}$ is not necessarily small, even with a tight spatio-temporal window. In such cases $s_{\mathrm{12}}^{\mathrm{2}}$ will be larger than a combined uncertainty 2σ², thus the estimate ${\widehat{\mathit{\sigma }}}^{\mathrm{2}}=s_{\mathrm{12}}^{\mathrm{2}}/\mathrm{2}$ will be biased high. Second, the number of self-collocated measurements for limb satellites is limited. Self-collocated measurements are usually found around the Poles, while at other latitudes a larger temporal separation (which can involve measurements in different day/night conditions) has to be accepted. In addition, the number of self-collocated measurements for instruments with coarse sampling (stellar and solar occultation) is relatively low. For example, ∼ 200 collocated occultations (with spatial separation Δr less than 300 km and temporal separation Δt less than 3h) per year of the star S30 (notation in the GOMOS catalogue) can be found for GOMOS; all located near the North Pole in winter. For ACE-FTS, fewer than 100 self-collocations per year with the criteria $\mathrm{\Delta }t=\mathrm{3}\mathrm{h},\mathrm{\Delta }r=\mathrm{300}\mathrm{km}$ are found, and ∼ 400 self-collocated measurements per year can be found with the collocation criteria $\mathrm{\Delta }t=\mathrm{5}\mathrm{h},\mathrm{\Delta }r=\mathrm{500}\mathrm{km}$. For SAGE II and SAGE III/ISS, there are no self-collocated measurements with the abovementioned collocation criteria.

Provided many collocated measurements from the same instrument are available (self-collocations), the precision of the dataset can also be estimated by computing a so-called structure function D(ρ) (e.g., Tatarskii, 1961), or the RMS difference of the field as a function of increasing separation in time and in space:

$$\begin{array}{}\text{(5)}& D\left(\mathit{\rho }\right)=D({\mathit{r}}_{\mathrm{1}}-{\mathit{r}}_{\mathrm{2}})={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}〈{\left[f\left({\mathit{r}}_{\mathrm{1}}\right)-f\left({\mathit{r}}_{\mathrm{2}}\right)\right]}^{\mathrm{2}}〉\end{array}$$

where r₁ and r₂ are two locations and a vector $\mathit{\rho }=r_{\mathrm{1}}-r_{\mathrm{2}}$ is their spatio-temporal separation. In geostatistics, D is called the variogram (Cressie, 1993; Matheron, 1963; Wackernagel, 2003). When using experimental (noisy) data for evaluation of the variogram/structure function, the difference of an atmospheric parameter in two locations is defined not only by the natural variability of this atmospheric parameter, but also by uncertainty of the measurements. Therefore, with the spatio-temporal separation ρ→0, D(ρ) tends toward the random uncertainty variance ${\mathit{\sigma }}_{\mathrm{noise}}^{\mathrm{2}}$ (the offset at zero is called “nugget” in geostatistics). Since self-collocated measurements are from the same instrument, no biases between them are expected. Figure 2 illustrates the structure function method, which is discussed in details in Sofieva et al. (2021) and applied to TROPOMI total ozone measurements.

Figure 2 Reproduced from (Sofieva et al., 2021): The schematic representation of the structure function estimated from noisy measurements, ρ denotes spatio-temporal separation. Blue line: the structure function for noise-free data, red line: the structure function for experimental (noisy) data.

Figure 3 Colored thin curves: experimental precision estimates $S_{\mathrm{12}}/\sqrt{\mathrm{2}}$ for different separation distances; thick curves with errorbars: the mean ex-ante uncertainty estimate and its standard deviation for propagated instrumental noise (magenta) and full random uncertainty (red); MIPAS self-collocations close to the North and South Poles in 2005–2011 during local summer are used.

For ozone profiles from limb instruments, the structure function method is difficult to apply, as it requires a substantial number of measurements with close separation. An analogous method – evaluation of the one-dimensional structure function in polar regions (with transformation of temporal mismatch to spatial separation using the ECMWF wind field) – has been applied for validation of random uncertainty estimates of the MIPAS and GOMOS ozone profiles (Laeng et al., 2015; Laeng and Von Clarmann, 2021; Sofieva et al., 2014).

Figure 3 illustrates the application of the structure function method to MIPAS v8 ozone profiles, which have a detailed error characterization. With decreasing separation distance between measurements, ex-post uncertainties $s_{\mathrm{12}}/\sqrt{\mathrm{2}}$ approach to a curve, which is between ex-ante estimates for total random error (thick red curves) and for propagated instrumental noise (thick magenta curves).

3.4 Using measurements from different instruments

3.4.1 Method of Fioletov et al. (2006)

Fioletov et al. (2006) have proposed estimating simultaneously the random data uncertainties and natural variability from sample variances of two perfectly collocated datasets and the variance of their difference. We reproduce the formulae here, as we assess the application of this method. The Fioletov method relies on sample variances $s_i^{\mathrm{2}}$ of the collocated data:

$$\begin{array}{}\text{(6)}& s_i^{\mathrm{2}}={\mathit{\sigma }}_{\mathrm{nat}}^{\mathrm{2}}+{\mathit{\sigma }}_i^{\mathrm{2}},i=\mathrm{1},\mathrm{2}\end{array}$$

and the variance of their difference (Eq. 4), which is reduced to:

$$\begin{array}{}\text{(7)}& s_{\mathrm{12}}^{\mathrm{2}}={\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}},\end{array}$$

by assuming ${\mathit{\sigma }}_{\mathrm{0},\mathrm{nat}}^{\mathrm{2}}\approx \mathrm{0}$. It is also assumed that the bias between datasets is the same for the selected sample.

In Eqs. (6) and (7), ${\mathit{\sigma }}_{\mathrm{nat}}^{\mathrm{2}}$ is natural variability and ${\mathit{\sigma }}_i^{\mathrm{2}}$ are measurement precisions. Solving Eqs. (6) and (7) for ${\mathit{\sigma }}_{\mathrm{nat}}^{\mathrm{2}}$, ${\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}$ and ${\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}$, we get their experimental estimates based on sample variance:

$$\begin{array}{}\text{(8)}& \begin{array}{l}{\widehat{\mathit{\sigma }}}_{\mathrm{nat}}^{\mathrm{2}}=\mathrm{0.5}\left(s_{\mathrm{1}}^{\mathrm{2}}+s_{\mathrm{2}}^{\mathrm{2}}-s_{\mathrm{12}}^{\mathrm{2}}\right)\\ {\widehat{\mathit{\sigma }}}_{\mathrm{1}}^{\mathrm{2}}=\mathrm{0.5}\left(s_{\mathrm{1}}^{\mathrm{2}}-s_{\mathrm{2}}^{\mathrm{2}}+s_{\mathrm{12}}^{\mathrm{2}}\right)\\ {\widehat{\mathit{\sigma }}}_{\mathrm{2}}^{\mathrm{2}}=\mathrm{0.5}\left(s_{\mathrm{2}}^{\mathrm{2}}-s_{\mathrm{1}}^{\mathrm{2}}+s_{\mathrm{12}}^{\mathrm{2}}\right)\end{array}\end{array}$$

The uncertainty of the natural variability and precision estimates given by Eq. (8) depend on uncertainty of sample variances, which depend, in turn, on sample variances themselves and the number of measurements. The estimates are thus only as accurate as the least accurate of these parameters. In approximation of large samples (when uncertainty of the sample variance can be approximated by Eq. 2), the variance of the estimates (8) can be expressed in terms of “true” natural variability and precision variances ${\mathit{\sigma }}_{\mathrm{nat}}^{\mathrm{2}}$, ${\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}$ and ${\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}$ as (using Eqs. 2, 6–8):

$$\begin{array}{}\text{(9)}& \begin{array}{rl}\mathrm{var}\left({\widehat{\mathit{\sigma }}}_{\mathrm{1}}^{\mathrm{2}}\right)& =\mathrm{var}\left({\widehat{\mathit{\sigma }}}_{\mathrm{2}}^{\mathrm{2}}\right)=\mathrm{var}\left({\widehat{\mathit{\sigma }}}_{\mathrm{nat}}^{\mathrm{2}}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}N}}\left({\left({\mathit{\sigma }}_{\mathrm{nat}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}\right)}^{\mathrm{2}}\right.\\ & \left.+{\left({\mathit{\sigma }}_{\mathrm{nat}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}\right)}^{\mathrm{2}}+{\left({\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}\right)}^{\mathrm{2}}\right)\end{array}\end{array}$$

with the following simple estimates for upper and lower limits (after opening brackets in Eq. 9):

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{\mathrm{1}}{N}}\left({\mathit{\sigma }}_{\mathrm{nat}}^{\mathrm{4}}+{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{4}}+{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{4}}\right)<\mathrm{var}\left({\widehat{\mathit{\sigma }}}_{\mathrm{1},\mathrm{2},\mathrm{nat}}^{\mathrm{2}}\right)<{\displaystyle \frac{\mathrm{1}}{N}}{\left({\mathit{\sigma }}_{\mathrm{nat}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}\right)}^{\mathrm{2}}.\end{array}$$

Since the precision estimates by the Fioletov method are linear combinations of three sample variances, they can have large uncertainty if one of the sample variances is large and/or the number of collocated measurements is limited. Not for all combinations of limb instruments perfectly collocated measurements can be found (especially for instruments with sparse sampling). In practice, satellite measurements separated by a few hundreds of kilometers and a few hours are considered collocated. The natural variability within the space-time collocation window is small but not zero, thus resulting in additional difficulties in the application of this method. Note that the estimates from Eq. (8) do not ensure positivity of ${\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}$ and ${\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}$. Negative solutions can be within uncertainty intervals; their appearance can be caused either by insufficient amount of data or by the unaccounted natural variability within the collocation window.

Figure 4 Application of Fioletov's method to MIPAS and SCIAMACHY ozone datasets in 2007. Left: sample standard deviations s₁ and s₂ in collocated pairs, and the standard deviation of differences s₁₂. Right: ex-post uncertainty estimates ${\widehat{\mathit{\sigma }}}_{\mathrm{1}}$ and ${\widehat{\mathit{\sigma }}}_{\mathrm{2}}$, and the estimate of natural variability ${\widehat{\mathit{\sigma }}}_{\mathrm{nat}}$ with 1σ uncertainties (solid lines). Ex-ante uncertainty estimates are shown with dashed lines. Error bars are 1σ uncertainties.

For illustration, we applied this method to MIPAS and SCIAMACHY measurements in 2007. Collocated profiles with time separation less than 5 h, spatial distance less than 400 km and latitude difference less than 2° were selected in the tropics (20° S–20° N), with 13 785 such profile pairs found. The left panel of Fig. 4 shows sample standard deviations s₁ and s₂ of MIPAS and SCIAMACHY profiles, respectively, and the standard deviation of differences s₁₂. The right panel shows the a posteriori (“ex-post” in terminology of von Clarmann et al., 2020) estimates of random uncertainties and natural variability from Eq. (8) with uncertainties therein given by Eq. (9). The estimates of random errors reported by the retrieval algorithms (“ex-ante” in terminology of von Clarmann et al., 2020) are also shown in right panels of Fig. 4 by dashed lines. Negative estimates of ${\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}$ and ${\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}$ are ignored. All computations are performed in absolute units, but the estimates are plotted as a percentage for clarity. We observe that the ex-ante and ex-post uncertainties of MIPAS ozone profiles are very close to each other. For SCIAMACHY, Fioletov's method suggests a larger uncertainty estimates at altitudes 25–37 km than reported in the retrievals.

Figure 5 Ex-ante and ex-post uncertainty estimates from Fioletov's method for MLS and ACE-FTS in years 2018–2019 (a) and 2018–2022 (b), and SAGE III/ISS and ACE-FTS in years 2018–2023. The number of collocated profiles is indicated in the panel titles.

As a general note for this and subsequent illustrations, the ex-ante random uncertainties for some instruments (see Table 1 for details) are due to measurement noise. This is a dominating source of random error, however not the only one (von Clarmann et al., 2020). As a result, ex-post random uncertainty estimates are expected to be slightly larger.

The best performance of the Fioletov's method is expected for datasets with dense sampling and similar random uncertainty estimates. The data should be selected in the regions of low variability. As mentioned above, Fioletov's method requires a large number of collocated profiles to yield reliable estimates of ex-post uncertainties. The application of this method to solar occultation data by ACE-FTS and SAGE III/ISS is illustrated in Fig. 5. The same collocation criteria are used, but the number of collocated profiles is significantly smaller than for MIPAS and SCIAMACHY, even though more years of data are used. For MLS and ACE-FTS, the number of collocations is 741 for years 2018–2019 and 1471 for years 2018–2022. The ex-post uncertainty estimates from Fioletov's method have substantial error bars (Fig. 5a, b). In application to ACE-FTS and SAGE III/ISS, only 19 collocated profiles are found, so the resulting uncertainty estimates have huge error bars (Fig. 5c).

The Fioletov's method is applicable to both individual-profile and tomographic retrievals, as the data in collocated pair can be considered as uncorrelated. As a general note, the fully rigorous approach of defining co-locations for tomographic retrievals should be different from that used for 1D individual-profile retrievals. Since the along-track dimension is part of the retrieval in the tomographic technique, interpolation should be performed along this dimension to the co-location point.

3.4.2 A differential method: comparisons of natural variability patterns

Sofieva et al. (2014) proposed a simple method for detecting flaws and/or checking of consistency of random uncertainty estimates. The authors called it a “differential method”. Let us consider, for example, two datasets selected in a region of small and slowly changing natural variability. A large sample size is assumed. If the random uncertainty estimates for both datasets are correct, then the difference in sample variance $s_{\mathrm{1}}^{\mathrm{2}}-s_{\mathrm{2}}^{\mathrm{2}}$ will be equal to the difference in precision estimates ${\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}-{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}$. The term ${\mathit{\sigma }}_{\mathrm{nat}}^{\mathrm{2}}$ from Eq. (3) cancels out because it is assumed to be the same for both samples. The estimates of the sample variance, $s_i^{\mathrm{2}}$, provide the upper limit for experimental estimates of measurement precision, as $s_i^{\mathrm{2}}>{\mathit{\sigma }}_i^{\mathrm{2}}$.

A simple comparison of sample variance $s_i^{\mathrm{2}}$ with the random uncertainty estimate ${\mathit{\sigma }}_i^{\mathrm{2}}$ enables the detection of overestimated random uncertainties, if the relation $s_i^{\mathrm{2}}>{\mathit{\sigma }}_i^{\mathrm{2}}$ is violated. Through such a comparison, Sofieva et al. (2014) found overestimated random uncertainties for the GOMOS ozone profiles using very dim stars, and further investigation by the instrument experts identified a flaw in accounting for instrumental dark charge noise.

If one of the datasets has realistic precision estimates, for example from well-calibrated instruments (so-called Fiducial Reference Measurements), or those estimates are validated by other methods, then application of the differential method is straightforward.

If there are several datasets with unvalidated (or not completely validated) uncertainty estimates, one can consider confronting natural variability estimates ${\widehat{\mathit{\sigma }}}_{\mathrm{nat}}^{\mathrm{2}}=s_i^{\mathrm{2}}-{\mathit{\sigma }}_i^{\mathrm{2}}$. Since the natural variability estimates from various datasets should agree within uncertainty intervals, strong deviations from the majority estimates can indicate potential flaws in the random error estimation.

For example, Sofieva et al. (2014) compared estimates of natural variability in the tropics from GOMOS data using different stars and found consistent positive values of ${\widehat{\mathit{\sigma }}}_{\mathrm{nat}}^{\mathrm{2}}$ for bright stars (${\mathit{\chi }}_{\mathrm{norm}}^{\mathrm{2}}\approx \mathrm{1}$ and application of the structure function method also suggested that the random uncertainties are realistic for bright stars). However, for very dim stars, negative values of ${\widehat{\mathit{\sigma }}}_{\mathrm{nat}}^{\mathrm{2}}$ have been detected, which, together with ${\mathit{\chi }}_{\mathrm{norm}}^{\mathrm{2}}<\mathrm{1}$, pointed to overestimated random uncertainties.

Figure 6 Top panels: Sample standard deviation s (solid lines) and the mean uncertainty estimates σ (dashed lines) in the tropical stratosphere (20° S–20° N) in 2002–2004 (1st column), 2006–2008 (2nd column) and 2018–2020 (3rd and 4th column). Bottom panels: the estimates of the natural variability ${\widehat{\mathit{\sigma }}}_{\mathrm{nat}}=\sqrt{s^{\mathrm{2}}-{\mathit{\sigma }}^{\mathrm{2}}}$ with its uncertainty (1 SD) in the tropics for the same periods. Colors are specified in the bottom panels. For GOMOS, occultations of 30 brightest stars are selected for the analysis.

In this paper, we illustrate this differential method by considering sample variance and uncertainty estimates from several limb and occultation instruments. The measurements are selected in the tropics, 20° S–20° N, in three periods, 2002–2004 (first column of Fig. 6), 2006–2008 (2nd column) and 2018–2020 (3rd and 4th column of Fig. 6), for the limb instruments operating in these periods. The upper panels of Fig. 6 show the sample standard deviation (solid lines) and the mean random uncertainty estimates (dashed lines). The lower panels show the estimates of natural variability ${\widehat{\mathit{\sigma }}}_{\mathrm{nat}}$ with associated uncertainties indicated by error bars. For GOMOS, occultations of the 30 brightest stars are used for the analysis, in order to make the GOMOS dataset more homogeneous and to avoid data with overestimated uncertainties. SAGE II and SAGE III/ISS ozone profiles were smoothed down to 2 km vertical resolution, for compatibility with other datasets.

As observed in Fig. 6, for all datasets except for OMPS USask, s<σ as expected. In case of OMPS USask, the mean uncertainty σ exceeds the sample standard deviation s at several altitudes, which indicates an overestimation of the random uncertainty component. The overestimation is caused by a bug in the V1.3.0 product and will be fixed in a future version. The profiles of natural variability ${\widehat{\mathit{\sigma }}}_{\mathrm{nat}}$ obtained from GOMOS, MIPAS, OSIRIS, ACE-FTS, SAGE II, SAGE III/ISS, OMPS UBr and OMPS NASA are very close to each other, and their uncertainty intervals overlap. For SCIAMACHY, the pattern of natural variability is also similar, but the increased sample variance at 25–35 km is not explained by its random uncertainty estimates; this suggests a slight underestimation of random error component at these altitudes. For MLS and SAGE III/ISS, there is a very good agreement with other datasets below 40–45 km; above 45 km, the random uncertainty estimates grow fast with altitude, which results in somewhat smaller estimates of ${\widehat{\mathit{\sigma }}}_{\mathrm{nat}}$ compared to other datasets. This indicates an overestimation of MLS and SAGE III/ISS random uncertainties above 40–45 km (see also the explanation in Sect. 2.4).

Successful application of this method implies the following conditions:

• a.

  Natural variability should be the same for both samples.

• b.

  Natural variability should not be too large compared to the precision estimates, otherwise the sample variance estimates will have large uncertainty. This condition of small natural variability is satisfied for ozone in the tropical stratosphere and in the summer stratosphere at other latitudes. (However, the random uncertainty estimates for limb-instruments are usually smaller than the natural variability.)

• c.

  Measurements within each sample should have similar precision.

The method can also be applied to the data retrieved with the tomographic approach, if the selected region is sufficiently large (exceeding the horizontal correlation length of tomographic retrievals).

If the natural variability is known from an external source (for example, estimated from the measurements with realistic uncertainty estimates or from a model that correctly reproduces natural variability), ex-post uncertainties can be estimated as ${\widehat{\mathit{\sigma }}}_{\text{ex-post}}^{\mathrm{2}}=s^{\mathrm{2}}-{\mathit{\sigma }}_{\mathrm{nat}}^{\mathrm{2}}$, where s² is the sample variance of a set of measurements and ${\mathit{\sigma }}_{\mathrm{nat}}^{\mathrm{2}}$ is the estimated of the natural variability. The use of the modelled data in validation of random uncertainty estimates is also discussed in Sect. 3.5 of this paper.

3.4.3 Triple collocation methods

Stoffelen's method

The idea of using the collocated measurements from three (or more) systems for data calibration and validation of uncertainties was proposed by Stoffelen (1998). In his formulation, it is supposed that three measurement systems X, Y, Z provide collocated measurements of the same quantity t. Let system X be the reference system with respect to which systems Y and Z are to be calibrated. Suppose also that linear calibration (simple scaling) is sufficient for the whole range of values under consideration, and that the reference system X is free of bias. Then the measurements can be written as

$$\begin{array}{}\text{(11)}& \begin{array}{l}x=t+{\mathit{\epsilon }}_x\\ y=c_y(t+{\mathit{\epsilon }}_y)\\ z=c_z(t+{\mathit{\epsilon }}_z),\end{array}\end{array}$$

where c_y and c_z are scaling factors and ${\mathit{\epsilon }}_x,{\mathit{\epsilon }}_y,{\mathit{\epsilon }}_z$ are random errors in each measurement sample. The random error components are assumed to be unbiased and not correlated with each other and with the parameter t. The calibration coefficients can be derived from covariances

$$\begin{array}{}\text{(12)}& \begin{array}{l}{c}_y=\mathrm{cov}(y,z)/\mathrm{cov}(x,z)\\ c_z=\mathrm{cov}(y,z)/\mathrm{cov}(x,y),\end{array}\end{array}$$

where cov$(\cdot ,\cdot )$ denotes covariance. These coefficients allow creating the calibrated data $y^{\ast }=c_y^{-\mathrm{1}}y,z^{\ast }=c_z^{-\mathrm{1}}z$. Then the natural variability of the parameter t can be estimated as e.g. ${\mathit{\sigma }}_t^{\mathrm{2}}=\mathrm{cov}(x,y^{\ast })$ and uncertainty variances as

$$\begin{array}{}\text{(13)}& \begin{array}{l}{\mathit{\sigma }}_x^{\mathrm{2}}=\mathrm{var}\left(x\right)-{\mathit{\sigma }}_t^{\mathrm{2}}\\ {\mathit{\sigma }}_y^{\mathrm{2}}=\mathrm{var}\left(y^{\ast }\right)-{\mathit{\sigma }}_t^{\mathrm{2}}\\ {\mathit{\sigma }}_z^{\mathrm{2}}=\mathrm{var}\left(z^{\ast }\right)-{\mathit{\sigma }}_t^{\mathrm{2}}.\end{array}\end{array}$$

The main assumptions of the method are: (a) the measurements are a linear function of the true signal with additive zero-mean random measurement noise; (b) measurement errors and true signal are stationary, and they are independent; (c) measurement errors are independent, and (d) the measurements are perfectly collocated, i.e. mismatch uncertainty is zero. The assumptions (b–d) are similar to other methods described above. The application of this method to validation of random uncertainties of tropospheric ozone from nadir instruments can be found in Hubert et al. (2021). In case of limb satellite observations, the requirement of triple collocation dramatically reduces the sample size (by approximately an order of magnitude); this results in larger uncertainties of the estimated parameters.

Another variant of the triple collocation method is described in the following subsection.

Von Clarmann's method

Von Clarmann proposed a method for random uncertainty validation, which uses 3 sets of measurements with pairwise collocations (Laeng and von Clarmann, 2021). This method takes into account the small-scale natural variability, which is estimated using a high-resolution chemistry-transport model data. Let us denote by $s_{ij}^{\mathrm{2}}$ the sample variance of differences in collocated datasets i and j, by ${\mathit{\nu }}_{ij}^{\mathrm{2}}$ – natural variability (mismatch) variance of the collocated datasets i and j, and by σ_i- ex-ante uncertainty estimates. Let us assume that the true random uncertainties are ${\mathit{\sigma }}_{\mathrm{true},i}^{\mathrm{2}}=c_i{\mathit{\sigma }}_i^{\mathrm{2}},i=\mathrm{1},\mathrm{2},\mathrm{3}$. Then the expressions for the sample variance in the collocated pairs results in the following system for determination of correction factors c_i:

$$\begin{array}{}\text{(14)}& \begin{array}{c}{c}_{\mathrm{1}}{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}+c_{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}+{\mathit{\nu }}_{\mathrm{12}}^{\mathrm{2}}=s_{\mathrm{12}}^{\mathrm{2}}\\ c_{\mathrm{1}}{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}+c_{\mathrm{3}}{\mathit{\sigma }}_{\mathrm{3}}^{\mathrm{2}}+{\mathit{\nu }}_{\mathrm{13}}^{\mathrm{2}}=s_{\mathrm{13}}^{\mathrm{2}}\\ c_{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}+c_{\mathrm{3}}{\mathit{\sigma }}_{\mathrm{3}}^{\mathrm{2}}+{\mathit{\nu }}_{\mathrm{23}}^{\mathrm{2}}=s_{\mathrm{23}}^{\mathrm{2}}\end{array}\end{array}$$

If ${\mathit{\nu }}_{ij}^{\mathrm{2}}$ are known, the solution of the linear system (14) is

$$\begin{array}{}\text{(15)}& \begin{array}{c}{c}_{\mathrm{1}}=\frac{\mathrm{1}}{\mathrm{2}{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}}\left[(s_{\mathrm{12}}^{\mathrm{2}}-{\mathit{\nu }}_{\mathrm{12}}^{\mathrm{2}})+(s_{\mathrm{13}}^{\mathrm{2}}-{\mathit{\nu }}_{\mathrm{13}}^{\mathrm{2}})-(s_{\mathrm{23}}^{\mathrm{2}}-{\mathit{\nu }}_{\mathrm{23}}^{\mathrm{2}})\right]\\ c_{\mathrm{2}}=\frac{\mathrm{1}}{\mathrm{2}{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}}\left[(s_{\mathrm{12}}^{\mathrm{2}}-{\mathit{\nu }}_{\mathrm{12}}^{\mathrm{2}})+(s_{\mathrm{23}}^{\mathrm{2}}-{\mathit{\nu }}_{\mathrm{23}}^{\mathrm{2}})-(s_{\mathrm{13}}^{\mathrm{2}}-{\mathit{\nu }}_{\mathrm{13}}^{\mathrm{2}})\right]\\ c_{\mathrm{3}}=\frac{\mathrm{1}}{\mathrm{2}{\mathit{\sigma }}_{\mathrm{3}}^{\mathrm{2}}}\left[(s_{\mathrm{13}}^{\mathrm{2}}-{\mathit{\nu }}_{\mathrm{13}}^{\mathrm{2}})+(s_{\mathrm{23}}^{\mathrm{2}}-{\mathit{\nu }}_{\mathrm{23}}^{\mathrm{2}})-(s_{\mathrm{12}}^{\mathrm{2}}-{\mathit{\nu }}_{\mathrm{12}}^{\mathrm{2}})\right]\end{array}\end{array}$$

The Eqs. (14) and (15) are written in terms of correction factors, as it is presented in the original report; however, they can be also presented in terms of ex-ante and ex-post uncertainties. Similarly to Fioletov's method, Eq. (15) does not guarantee positivity of c_i (solution of the linear system), which may result in unphysical negative estimates of random error variance. Since the estimates by Eq. (15) are the linear combinations of sample variances, the accurate estimate require large samples of collocated data, similarly to the Fioletov method.

For illustration, the method was applied to MIPAS, MLS, and SCIAMACHY data in 2007 at 20° S–20° N, where 3236 triple collocations (with time difference <4 h and spatial separation <300 km) are found. The small-scale variability estimates were obtained from BASCOE model data field down-sampled to typical horizontal resolution along line of sight of limb instrument data (Laeng et al., 2022).

The estimation of the small-scale variability using the ozone fields generated by the advanced chemistry-transport model with a high horizontal resolution seems to be correct, as this mismatched variability is caused by dynamics, and it is characterized by the statistical characteristic (variance). Although some small-scale processes may not be fully resolved in the model, this should not influence the ν_ij estimates, as the effective horizontal resolution of the limb measurements along the line of sight is ∼ 300 km.

The results of application of the von Clarmann's method (Fig. 7) agree very well with other methods presented in our paper (see also below).

Figure 7 Dashed lines: ex-ante uncertainties; solid lines with error bars: ex-post random uncertainties estimated by the von Clarmann method. Error bars are 1σ uncertainties.

3.5 Using CTM simulation in validation of uncertainties

Modern chemistry-transport models have high horizontal and vertical resolution, and a majority of them use meteorological reanalyses in the advection schemes. They show good agreement with the observational data, therefore it is attractive to use the information the models provide in validation of uncertainties. For example, the modelled field can be used for characterization of differences due to the co-location mismatch, i.e. differences in spatio-temporal sampling and smoothing of the variable and inhomogeneous ozone field. Such an approach has been applied in several studies (e.g., Sheese et al., 2021; Verhoelst et al., 2015).

The model estimates of small-scale natural variability are used also in von Clarmann's method. Potentially, analogous characterization would also improve the Fioletov's method.

Sofieva et al. (2022) used ozone data, which are simulated with the chemistry-transport model SILAM adjusted to MLS for ex-post random uncertainty estimates by the differential method (this method is referred to as the SUNLIT method). For each instrument and each month, the authors evaluated sample variance in 10° latitude zones from experimental data and the SILAM-adjusted field, which is sub-sampled at measurements locations. The sample variance of the model dataprovides the estimates of natural variability. Then ex-post uncertainties are estimated as ${\widehat{\mathit{\sigma }}}_{\text{ex-post}}^{\mathrm{2}}=s^{\mathrm{2}}-{\mathit{\sigma }}_{\mathrm{nat}}^{\mathrm{2}}$, where s² is the sample variance in a set of measurements and ${\mathit{\sigma }}_{\mathrm{nat}}^{\mathrm{2}}$ is the estimate of the natural variability. Figure 8 illustrates ex-ante and ex-post uncertainties for GOMOS, MIPAS, SCIAMACHY, OSIRIS and MLS using the data in September 2007. These estimates agree well with those obtained by the Fioletov's and von Clarmann's method (Fig. 4, right, Figs. 7, and 9). The approach of Sofieva et al. (2022) allows selecting sufficiently large data samples in a relatively short time period; the authors applied their method to the adjustment of random uncertainties for each month and for each instrument.

Figure 8 Ex-ante (dashed lines) and ex-post (solid lines) random uncertainty estimates for GOMOS, MIPAS, SCIAMACHY, OSIRIS and MLS in September 2007, based on the method described in Sofieva et al. (2022).

Figure 9 Uncertainty estimates for MIPAS at 20° S–20° N in 2007 by different methods. Dashed lines: ex-ante uncertainties (coinciding for Fioletov's and von Clarmann's methods), solid lines: ex-post uncertainties (%). For the SUNLIT method, data in September 2007 are used. Error bars are 1σ uncertainties.

Figure 9 compares the ex-post uncertainties for MIPAS ozone retrievals in 2007 at latitudes 20° S–20° N evaluated by different methods. For Fioletov's and von Clarmann's methods, the data from the entire year 2007 are used, while for SUNLIT only data from September 2007 are used. The agreement between ex-post uncertainty estimates by different methods is good: ex-post uncertainties are within the corresponding error bars. In the upper stratosphere, the SUNLIT method yields slightly larger ex-post uncertainties, probably due to underestimated natural variability in the model data (in SUNLIT, daily mean ozone profiles are used, thus missing diurnal variability).

3.6 Notes on validation with fiducial reference measurements

If a dataset with well known (or validated and concluded to be realistic) uncertainty estimates is available, then the validation of uncertainties of a second dataset using either Eq. (4) or the differential method is straightforward. Such an approach is usually used in validation of satellite measurements with ground-based data. When exploiting Eq. (4), it is advantageous to characterize/simulate ${\mathit{\sigma }}_{\mathrm{0},\mathrm{nat}}^{\mathrm{2}}$. Such approach was explored for the validation of satellite total ozone column data by ground-based measurements (Verhoelst et al., 2015).

For ozone profiles, ozonesondes are usually used for validation of satellite data (evaluation of biases and drifts). However, according to Tarasick et al. (2021) the characterization of ozonesonde uncertainties is even more complicated than for satellite data, and the uncertainties are not constant but varying in the range 5 %–10 % (sometimes even up to 20 %). Together with a limited number of tight collocations, these features impose limitations for validation of satellite random uncertainties using ozonesonde data.

3.7 On using the Markov Chain Monte Carlo method

The random uncertainties reported by retrieval algorithms are usually estimated via propagation of instrumental noise and other random uncertainties through the inversion algorithm. Markov chain Monte Carlo (MCMC) method can be used to produce a robust estimate of the probability distribution of a retrieved quantity that is nonlinearly related to the measurements and that has non-Gaussian error statistics. A methodology for validating the traditional error characterization by applying the MCMC technique can be found in e.g. Tamminen (2004). This paper shows the application of MCMC method to GOMOS data. The MCMC technique is suitable for studying uncertainties of retrieved parameters and it enables analyzing the error structure also in a nonlinear case (and thus validating the standard Gaussian characterization). The advantage of the sampling based MCMC method is also that it allows implementing non-Gaussian measurement and modelling error characterization as well as using non-Gaussian prior information. While the MCMC method cannot provide information about missing or overestimated uncertainties directly, the method is often implemented so that unknown uncertainties are parametrized and included, e.g., via hierarchical formulation, allowing these uncertainties to be taken into account. It would be very useful to compare such approaches with the ones presented in this overview paper, in the future.

4 Summary and discussion

In this paper, we presented methods for random uncertainty validation and illustrated their application using ozone profiles retrieved from measurements by satellite instruments in the limb-viewing geometry. These methods considered in this study rely on deriving a posteriori (ex-post) random uncertainties using statistical analyses of collocated data samples. Advantages and limitations of each method are discussed, as well as accuracy of ex-post random uncertainty estimates.

As a general requirement for all methods, the data samples should be selected in regions of small and slowly changing natural variability. Otherwise, if the natural variability exceeds significantly random uncertainties, this prevents the computation of reliable ex-post estimates of random uncertainties. The methods for random uncertainty validation are divided into two groups depending on what kind of data are used: (1) from the same instrument and (2) from different instruments.

Practical examples of validation of random uncertainty with the discussion of advantages and limitations of each method are provided in this study. It is shown that, for instruments with dense sampling, such as MIPAS and MLS, several methods can be applied, for example those based on self-collocations or collocations with other datasets. For datasets that are obtained with tomographic retrievals, the Fioletov's, von Clarmann's and differential methods can be applied. For instruments with coarse sampling, such as GOMOS, ACE-FTS or SAGE II-III, the differential method is the most appropriate. It has been shown previously and also confirmed in this study that simulations with high-quality and high-resolution chemistry-transport models are useful in validation of reported random uncertainties: the model simulations can be used for estimation of small-scale natural variability.

The methods presented in this overview can be also applied to other measurements. In particular, the structure function method has been already successfully applied to total column measurements by TROPOMI in Sofieva et al. (2021). All methods can be applied also to data with coarse vertical resolution, such as profiles retrieved from nadir-looking instruments. For the application of the methods based on the statistics of differences, the profiles should have a compatible vertical resolution. This might require prior application of harmonization (see Keppens et al., 2019 for details). Then the validation of random uncertainties can be performed at the vertical scales corresponding to harmonized profiles.