The IASI NH₃ version 4 product: averaging kernels and improved consistency

Satellite measurements play an increasingly important role in the study of atmospheric ammonia (NH₃). Here, we present version 4 of the Artificial Neural Network for IASI (ANNI; IASI: Infrared Atmospheric Sounding Interferometer) retrieval of NH₃. The main change is the introduction of total column averaging kernels (AVKs), which can be used to undo the effect of the vertical profile shape assumption of the retrieval. While the main equations can be matched term for term with analogous ones used in UV/Vis retrievals for other minor absorbers, we derive the formalism from the ground up, as its applicability to thermal infrared measurements is non-trivial. A large number of other smaller changes were introduced in ANNI v4, most of which improve the consistency of the measurements across time and across the series of IASI instruments. This includes a more robust way of calculating the hyperspectral range index (HRI), explicitly accounting for long-term changes in CO₂ in the HRI calculation and the use of a reprocessed cloud product that was specifically developed for climate applications. The NH₃ distributions derived with ANNI v4 are very similar to the ones derived with v3, although values are about 10 %–20 % larger due to the improved setup of the HRI. We exclude further large biases of the same nature by showing the consistency between ANNI v4 derived NH₃ columns with columns obtained with an optimal estimation approach. Finally, with v4, we revised the uncertainty budget and now report systematic uncertainty estimates alongside random uncertainties, allowing realistic mean uncertainties to be estimated.

1 Introduction

Atmospheric ammonia (NH₃) primarily originates from agriculture and related activities. Its presence in the atmosphere leads to a reduction of life quality and to millions of premature deaths via its contribution to particulate matter (Pozzer et al., 2017). As one of the main forms of reactive nitrogen (Galloway et al., 2021), NH₃ is also a key element in the global nitrogen cycle with devastating effects on the environment when deposited in excess (Sutton et al., 2014).

Satellite measurements of NH₃ abundances have in the past decade contributed to our understanding of its global distribution, spatiotemporal variations, emission sources, concentration trends, transport patterns, chemistry and deposition levels. Currently, the two most widely used satellite datasets are those derived from observations of the Cross-track Infrared Sounder (CrIS) (Shephard et al., 2020) and the three Infrared Atmospheric Sounding Interferometers (IASI) (Van Damme et al., 2021). The CrIS product relies on optimal estimation, while the IASI product is based on the conversion of a spectral NH₃ index to a total column.

The first version of the IASI-NH₃ product (Van Damme et al., 2014) used look-up-tables (LUT) for the conversion. The LUTs were replaced with a more flexible neural network (NN) in Whitburn et al. (2016). Since then, the Artificial Neural Network for IASI (ANNI) retrieval approach underwent a series of incremental improvements that are documented in Van Damme et al. (2017), Franco et al. (2018) and Van Damme et al. (2021). In Franco et al. (2018, 2019, 2020, 2022) and Rosanka et al. (2021) the ANNI retrieval framework was expanded for the retrieval of other minor trace gases from IASI observations. A similar retrieval approach was also recently adopted for isoprene retrievals from CrIS (Wells et al., 2022).

Returning to ANNI-NH₃, in Van Damme et al. (2017), a reanalysis product was introduced. This product differs in the origin of the input parameters that are used for the NN. Whereas the baseline product (also called near-real-time (NRT) product) uses operational IASI Level 2 (L2) information on the pressure, humidity and temperature profiles, and a climatology characterizing the boundary layer height, the reanalysed product uses ERA5 model output for these parameters (Hersbach et al., 2020), interpolated in time and space to match the observations. The resulting product is temporally more consistent as it removes the effect of the several changes that occurred in the L2 products throughout the years. Both NH₃ products include several empirical corrections to counter small differences observed between the three IASI instruments and small biases that occurred as a result of sporadic changes to the IASI instrument or in the L0 to L1c processing of the spectra.

The IASI NH₃ product is widely used by the scientific community. However, the absence of averaging kernels (AVKs) has in the past hampered model comparison and assimilation. This paper presents version 4 of the ANNI retrieval framework, with the most important change being the introduction of AVKs. After a brief recapitulation of the retrieval algorithm in Sect. 2, the AVK framework is presented in Sect. 3. This includes its theoretical basis, practicalities related to how the AVKs are calculated within ANNI and a discussion on how they can be used in measurement-model comparison and assimilation. Other changes that were introduced in ANNI v4 are detailed in Sect. 4 (for those related to temporal consistency) and Sect. 5 (for all other changes). In Sect. 6, an evaluation of the NH₃ product is presented, comparing the v3 with the v4 product and the neural network output with retrievals based on optimal estimation. In the final part of this paper, we present the revised uncertainty budget of the ANNI retrieval.

2 ANNI retrieval overview

Here, we give a brief overview of the ANNI algorithm and refer to the previously cited papers on the NH₃ product for a detailed description and the rationale behind the different retrieval choices. Table 1 summarizes the most important quantities and symbols used in this paper.

Table 1 Summary of the main quantities and associated symbols that are used in the AVK formalism.

The retrieval consists of two independent computational steps. The first one characterizes the NH₃ signal strength in a spectrum L, via the so-called hyperspectral range index (HRI), which relies on a mean spectrum $\stackrel{\mathrm{‾}}{\mathit{L}}$ and associated covariance matrix S constructed from a set of spectra with no observable NH₃ spectral signatures. It is defined as

$$\begin{array}{}\text{(1)}& \text{HRI}=N{\mathit{K}}^{\mathrm{T}}{\mathbf{S}}^{-\mathrm{1}}\left(\mathit{L}-\stackrel{\mathrm{‾}}{\mathit{L}}\right),\end{array}$$

with N a normalization constant and K an NH₃ Jacobian. By construction, the HRI has a mean of zero on the spectra from which $\stackrel{\mathrm{‾}}{\mathit{L}}$ and S are constructed. The normalization factor guarantees that the HRI has a standard deviation of 1 spectra containing only background levels of NH₃.

The second part of the algorithm relies on a neural network to link the measured HRIs to estimates ${\widehat{X}}^a$ of the true NH₃ total columns X, via a scaling factor SF^a:

$$\begin{array}{}\text{(2)}& {\widehat{X}}^a={\displaystyle \frac{\text{HRI}}{{\text{SF}}^a}}+B,\end{array}$$

with B an NH₃ background column corresponding to $\stackrel{\mathrm{‾}}{\mathit{L}}$. As in previous versions of ANNI, we will assume a zero background column for NH₃ for ANNI v4. However, in what follows, we develop the theory for an arbitrary background column B, so that the recipe can be applied to the other tracers retrieved from IASI with ANNI (e.g. CH₃OH or PAN). Retrieved quantities will be indicated with a hat, as in $\widehat{X}$. Superscripts refer to the assumed or modelled vertical profile shape. In the ANNI retrieval framework, the scaling factor SF^a is the quantity that is calculated with an NN, for each individual observation, based on the state of the atmosphere (temperature and water vapour profile, surface pressure), the surface temperature and emissivity, the zenith angle, the HRI and an assumed vertical profile shape. For NH₃, the volume mixing ratio (VMR) vertical profile is modelled as a Gaussian,

$$\begin{array}{}\text{(3)}& \text{VMR}\left(z\right)=\text{VMR}\left(z_{\mathrm{0}}\right)e^{-\frac{(z-z_{\mathrm{0}}{)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}},\end{array}$$

with z the altitude about ground level, z₀ the peak altitude and σ the width of the profile. Over land, the peak altitude is set at the surface, with a width σ equal to the boundary layer height. Over ocean, the peak altitude is set to 1.4 km with a σ of 0.9 km. In general, ${\widehat{X}}^a\ne X$ because of instrumental noise, errors in the assumed vertical profile, imperfect knowledge of one of the other input parameters and errors in the spectroscopic parameters or forward model. The NN is trained from a large set of forward modelled spectra. Appropriate pre- and post-retrieval flags accompany the retrieval. The pre-filter removes respectively measurements with erroneous L1 or excess cloud coverage. The post-filter flags retrievals with limited or no sensitivity to the measured quantity, satisfying

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\mathrm{1}}{\left|\text{SF}\right|}}>\mathrm{1.5}\times {\mathrm{10}}^{\mathrm{16}}\mathrm{molec}.{\mathrm{cm}}^{-\mathrm{2}},\end{array}$$

and retrievals whose HRIs are either too noisy or for which the assumed vertical profile is incompatible with the measured HRI, satisfying

$$\begin{array}{}\text{(5)}& \left|\text{HRI}\right|>\mathrm{1.5}\text{and}{\widehat{X}}^a<\mathrm{0}.\end{array}$$

Note that this filter only removes a fraction of the negative columns (see Clarisse et al., 2019; Whitburn et al., 2016 for a discussion on why it is important to keep these). Via propagation of uncertainty, a total retrieval uncertainty is calculated for each individual measurement alongside the retrieved column.

3 Averaging kernels

The general AVK framework that we introduce below bears a lot of similarity to the total column AVK formalism (Eskes and Boersma, 2003) developed for the DOAS retrieval approach of weakly absorbing species (see also Palmer et al., 2001; Rodgers and Connor, 2003; Boersma et al., 2004, 2016; Cooper et al., 2020). In the DOAS retrieval approach, the total column X is retrieved as

$$\begin{array}{}\text{(6)}& {\widehat{X}}^a={\displaystyle \frac{\text{SCD}}{{\text{AMF}}^a}},\end{array}$$

with SCD the slant column density and AMF the air mass factor which accounts for the atmospheric conditions and assumed vertical profile. Equation (6) has the same functional form as the main formula (Eq. 2) of the ANNI retrieval formalism, with the SCD corresponding to the HRI and the AMF to the SF provided by an NN.

One key element on which the total column AVK formalism of Eskes and Boersma (2003) relies is linearity and additivity of the spectrum with respect to changes in the trace gas amount. Linearity is a consequence of the curve of growth of spectral lines for low optical depths (Thorne et al., 1999) and in the DOAS approach also implies that the SCD is proportional to the trace gas abundance. Additivity represents the fact that the effect of different atmospheric layers can be summed up independently from each other. Given the definition of the HRI and the effects of thermal emission of the atmosphere, it is not obvious that these hold in the infrared spectral domain, and for this reason, we derive both properties below.

3.1 On the linearity and additivity of the HRI

Let $L_{\mathit{\nu }}^B$ be the radiance at the sensor for a scene with climatological background levels B of the target trace gas and a corresponding HRI of zero. Dividing the atmosphere in n appropriately spaced layers z, we denote by B_z the corresponding partial columns ($B={\sum }_z{B}_z$). We can then calculate $L_{\mathit{\nu }}^B$ from the following sequence of equations: (Rodgers, 2000; Petty, 2006)

$$\begin{array}{}\text{(7)}& {\displaystyle }{L}_{\mathrm{1}\mathit{\nu }}^B=t_{\mathrm{1}\mathit{\nu }}(L_{\mathrm{0}\mathit{\nu }}-B_{\mathrm{1}\mathit{\nu }})+B_{\mathrm{1}\mathit{\nu }}\text{(8)}& {\displaystyle }{L}_{\mathrm{2}\mathit{\nu }}^B=t_{\mathrm{2}\mathit{\nu }}(L_{\mathrm{1}\mathit{\nu }}-B_{\mathrm{2}\mathit{\nu }})+B_{\mathrm{2}\mathit{\nu }}\text{(9)}& {\displaystyle }\begin{array}{rl}& \mathrm{⋮}\\ & L_{\mathit{\nu }}^B=t_{n\mathit{\nu }}\left(L_{(n-\mathrm{1})\mathit{\nu }}-B_{n\mathit{\nu }}\right)+B_{n\mathit{\nu }}.\end{array}\end{array}$$

Forward substitution of the above relations yields $L_{\mathit{\nu }}^B$ as a function of the surface term L_0ν. Here t_zν are the layer transmittances that account for the absorption of all atmospheric species. The z dependence of this parameter is related to vertical variations in the atmospheric constituents and the pressure and temperature dependence of the line intensities. B_zν=B_ν(T_z) corresponds to the Planck's blackbody function for an averaged layer temperature T_z.

Figure 1 Numerical demonstration of the linearity and additivity of the HRI as a function of a change in partial column. In the blue and green scenario, NH₃ was varied at one fixed altitude. In the red scenario, partial columns in both layers were varied simultaneously. The solid black lines represent linearity, whereas the dash-dotted line, being the sum of the green line and blue line, represents additivity.

For this scene, we now introduce an additional trace amount X−B and write the corresponding observed radiance as L_ν. In each layer, the transmittance will decrease by a factor $t_{z\mathit{\nu }}^{X_z-B_z}=e^{-{\mathit{\tau }}_{z\mathit{\nu }}}\approx \mathrm{1}-{\mathit{\tau }}_{z\mathit{\nu }}$, with τ_zν the small optical depth caused by the excess X_z−B_z. With this, the sequence of equations becomes

$$\begin{array}{}\text{(10)}& {\displaystyle }{L}_{\mathrm{1}\mathit{\nu }}=t_{\mathrm{1}\mathit{\nu }}\left(\mathrm{1}-{\mathit{\tau }}_{\mathrm{1}\mathit{\nu }}\right)\left(L_{\mathrm{0}\mathit{\nu }}-B_{\mathrm{1}\mathit{\nu }}\right)+B_{\mathrm{1}\mathit{\nu }}\text{(11)}& {\displaystyle }{L}_{\mathrm{2}\mathit{\nu }}=t_{\mathrm{2}\mathit{\nu }}\left(\mathrm{1}-{\mathit{\tau }}_{\mathrm{2}\mathit{\nu }}\right)\left(L_{\mathrm{1}\mathit{\nu }}-B_{\mathrm{2}\mathit{\nu }}\right)+B_{\mathrm{2}\mathit{\nu }}\text{(12)}& {\displaystyle }\begin{array}{rl}& \mathrm{⋮}\\ & L_{\mathit{\nu }}=t_{n\mathit{\nu }}\left(\mathrm{1}-{\mathit{\tau }}_{n\mathit{\nu }}\right)\left(L_{(n-\mathrm{1})\mathit{\nu }}-B_{n\mathit{\nu }}\right)+B_{n\mathit{\nu }}.\end{array}\end{array}$$

For optical thicknesses of the target trace gas τ_zν well below one, the second (τ_iντ_jν) and higher order terms in optical depth can be neglected. Combining both sets of equations, one can verify that

$$\begin{array}{}\text{(13)}& {\displaystyle }{L}_{\mathrm{1}\mathit{\nu }}=L_{\mathrm{1}\mathit{\nu }}^B-{\mathit{\tau }}_{\mathrm{1}\mathit{\nu }}\left(L_{\mathrm{0}\mathit{\nu }}-B_{\mathrm{1}\mathit{\nu }}\right)t_{\mathrm{1}\mathit{\nu }}\text{(14)}& {\displaystyle }{L}_{\mathrm{2}\mathit{\nu }}=L_{\mathrm{2}\mathit{\nu }}^B-{\mathit{\tau }}_{\mathrm{1}\mathit{\nu }}\left(L_{\mathrm{0}\mathit{\nu }}-B_{\mathrm{1}\mathit{\nu }}\right)t_{\mathrm{1}\mathit{\nu }}{t}_{\mathrm{2}\mathit{\nu }}-{\mathit{\tau }}_{\mathrm{2}\mathit{\nu }}\left(L_{\mathrm{1}\mathit{\nu }}-B_{\mathrm{2}\mathit{\nu }}\right)t_{\mathrm{2}\mathit{\nu }}\text{(15)}& {\displaystyle }\begin{array}{rl}& \mathrm{⋮}\\ & L_{\mathit{\nu }}=L_{\mathit{\nu }}^B-\sum _{z=\mathrm{1}}^n{\mathit{\tau }}_{z\mathit{\nu }}\left(L_{(z-\mathrm{1})\mathit{\nu }}^B-B_{z\mathit{\nu }}\right)\prod _{j=z}^n{t}_{j\mathit{\nu }}.\end{array}\end{array}$$

Each of the terms in the sum expresses the effect of the absorption in one layer due to the excess X_z−B_z and the attenuation by the layers above. Note that the absorptions are proportional to both the optical depth and the local thermal contrast $\left(L_{(z-\mathrm{1})\mathit{\nu }}^B-B_{z\mathit{\nu }}\right)$, which are two parameters that drive the measurement sensitivity in the infrared (Bauduin et al., 2017). The optical thickness in turn is proportional to the partial column of the target species ${\mathit{\tau }}_{z\mathit{\nu }}\propto (X_z-B_z)$ and thus

$$\begin{array}{}\text{(16)}& L_{\mathit{\nu }}=L_{\mathit{\nu }}^B-\sum _z{c}_{z\mathit{\nu }}(X_z-B_z).\end{array}$$

The constants c_zν depend on the state of the atmosphere at level z and above, but are independent of the excess trace gas amount X−B.

The HRI is by definition a linear combination of spectral channels ν, from Eq. (1),

$$\begin{array}{}\text{(17)}& \text{HRI}=\sum _{\mathit{\nu }}{w}_{\mathit{\nu }}{L}_{\mathit{\nu }}+C,\end{array}$$

with C and w_ν numerical constants. Using Eq. (16) and the fact that the HRI is zero on the background ($\text{HRI}={\sum }_{\mathit{\nu }}{w}_{\mathit{\nu }}{L}_{\mathit{\nu }}^B+C=\mathrm{0}$), we find

$$\begin{array}{}\text{(18)}& {\displaystyle }\text{HRI}=\sum _{\mathit{\nu }}{w}_{\mathit{\nu }}{L}_{\mathit{\nu }}^B-\sum _{\mathit{\nu },z}{w}_{\mathit{\nu }}{c}_{z\mathit{\nu }}(X_z-B_z)+C\text{(19)}& {\displaystyle }=-\sum _{\mathit{\nu },z}{w}_{\mathit{\nu }}{c}_{z\mathit{\nu }}(X_z-B_z).\end{array}$$

This equation implies that for optically thin absorbers, the HRI can be written as a weighted sum of the partial column enhancements. This linearity and additivity of the NH₃ HRI as a function of partial columns is illustrated in Fig. 1 by means of simulations with a radiative transfer code. The blue line illustrates the nearly linear increase in HRI for increasing NH₃ columns in the first atmospheric layer of the model. Starting from about 7×10¹⁶ molec. cm⁻², slow departure from linearity is observed. The green line shows the same for the second atmospheric layer. The overall larger HRI in this second layer results from higher thermal contrast (TC) higher up in the atmosphere. Finally, the red line represents the HRI when NH₃ is increased simultaneously in both layers. The dash-dotted line is the sum of the blue and green curve and for low columns is almost indiscernible from the simulated HRI, illustrating additivity in the optically thin limit.

To make the link with Eq. (2), Eq. (19) can be rewritten as

$$\begin{array}{}\text{(20)}& {\displaystyle }\text{HRI}=\sum _z{\text{SF}}_z(X_z-B_z)\text{(21)}& {\displaystyle }=\sum _z{\text{SF}}_z{\displaystyle \frac{X_z-B_z}{X-B}}(X-B)\text{(22)}& {\displaystyle }=\text{SF}(X-B),\end{array}$$

with SF_z local scaling factors and

$$\begin{array}{}\text{(23)}& \text{SF}=\sum _z{\text{SF}}_z{\displaystyle \frac{X_z-B_z}{X-B}},\end{array}$$

total scaling factors that depend on the scene conditions (e.g. surface temperature, atmospheric temperature and pressure profiles, vertical profile of the target species). We note that whereas SF depends on the normalized vertical profile shape $\frac{X_z-B_z}{X-B}$, both SF_z and SF are independent of the total column for a fixed profile shape.

Finally, introducing

$$\begin{array}{}\text{(24)}& {\text{HRI}}_z={\text{SF}}_z(X_z-B_z),\end{array}$$

we see from Eq. (20) that the HRI can be decomposed in different partial HRI_z:

$$\begin{array}{}\text{(25)}& \text{HRI}=\sum _z{\text{SF}}_z(X_z-B_z)=\sum _z{\text{HRI}}_z.\end{array}$$

These HRI_z quantify how much each layer contributes to the total HRI and they are, again in the optically thin limit, independent from each other.

3.2 Total column averaging kernels

Equation (22) motivates the NN retrieval formula,

$$\begin{array}{}\text{(26)}& {\widehat{X}}^a={\displaystyle \frac{\text{HRI}}{{\text{SF}}^a}}+B,\end{array}$$

with SF^a an estimated scaling factor which depends on the best estimates of all the dependencies of SF, including the vertical profile of the target species. As we have just shown, in the optically thin limit SF^a is independent of the total column amount. This is no longer the case for strong absorptions when non-linear effects become increasingly important. In the ANNI retrieval framework this is taken care of by including the HRI as an input parameter in the calculation of the SF. In what follows, we will first derive the AVK formalism in the optically thin limit, assuming both linearity and additivity (and therefore SF that are independent of the (partial) column). In Sect. 3.5 we will then show how to correct for small errors that arise when these conditions are not met.

Figure 2 A priori vertical profile a_z, averaging kernel $A_z^a$ and vertical partitioning of signal $V_z^a$.

The quantities ${\widehat{X}}_z-B_z$ express the local enhancement of the trace gas column at an altitude z, and it is their relative proportions that are constrained in the retrieval by the use of an a priori profile shape. We write

$$\begin{array}{}\text{(27)}& a_z={\displaystyle \frac{{\widehat{X}}_z^a-B_z}{{\widehat{X}}^a-B}}\end{array}$$

for the normalized a priori vertical partial column profile, or partial column profile shape, with ${\widehat{X}}_z^a$ the partial columns corresponding to the retrieved ${\widehat{X}}^a$ and ${\sum }_z{a}_z=\mathrm{1}$. In the previous section we demonstrated that in the optically thin limit, the scaling factor SF^a equals the weighted sum of the different local scaling factors defined by the profile shape:

$$\begin{array}{}\text{(28)}& {\text{SF}}^a=\sum _z{\text{SF}}_z{\displaystyle \frac{X_z^a-B_z}{X^a-B}}=\sum _z{\text{SF}}_z{a}_z.\end{array}$$

Defining the averaging kernel as

$$\begin{array}{}\text{(29)}& A_z^a={\displaystyle \frac{{\text{SF}}_z}{{\text{SF}}^a}}={\displaystyle \frac{{\text{SF}}_z}{{\sum }_z{\text{SF}}_z{a}_z}},\end{array}$$

we can express ${\widehat{X}}^a$ as a function of X_z by eliminating the HRI from Eqs. (20) and (26):

$$\begin{array}{}\text{(30)}& {\displaystyle }{\widehat{X}}^a={\displaystyle \frac{\text{HRI}}{{\text{SF}}^a}}+B\text{(31)}& {\displaystyle }={\displaystyle \frac{{\sum }_z{\text{SF}}_z(X_z-B_z)}{{\text{SF}}^a}}+B\text{(32)}& {\displaystyle }=\sum _z{A}_z^a(X_z-B_z)+B.\end{array}$$

As can be seen from this equation, the averaging kernel $A_z^a$ fully characterizes the measurement, and can be used to mathematically map the true profile X_z to the measured total column ${\widehat{X}}^a$.

Two example AVKs for an NH₃ retrieval are shown in the middle panel of Fig. 2, for a retrieval over land and over ocean, with a priori profiles shown in the left panel. Both AVKs logically increase with altitude, as the temperature difference between the surface and a given atmospheric layer, and therefore the scaling factor, increases with altitude. Apart from a multiplicative constant $(\mathrm{1}/{\text{SF}}^a)$, the AVKs are independent from the a priori profiles, which explains why both land and ocean AVKs increase similarly with altitude. The multiplicative constant determines the altitude for which the AVK is one. This altitude can be interpreted as an equivalent effective NH₃ altitude, where all NH₃ can be thought to be located. For the ocean and land retrieval, this altitude is located respectively around 1.6 and 0.7 km, consistent with the a priori profiles shown in the left panel.

3.3 Interpretation

It is instructive to compare the total column averaging kernel, as defined above, with the one arising in optimal estimation retrievals (Rodgers, 2000):

$$\begin{array}{}\text{(33)}& {\widehat{\mathit{X}}}^{\mathit{a}}=\left(\mathbf{I}-\mathbf{A}\right){\mathit{X}}^{\mathit{a}}+\mathbf{A}\mathit{X}.\end{array}$$

Here, we use matrix and vector notation, with X and $\widehat{\mathit{X}}$ respectively the true and retrieved partial columns and A the AVK matrix. I is the identity matrix. Both Eqs. (32) and (33) allow simulating the retrieval process for any hypothetical X (e.g. from a model, or an independent measurement). However, this is largely where their similarity ends, as there are important differences when it comes to interpreting these two types of AVK.

The first difference is the role of the a priori. For the total column retrieval, the a priori fixes only the vertical profile shape, while for the optimal estimation retrieval, the a priori affects both the vertical profile shape and the retrieved value at each altitude separately. Equation (33) expresses that the retrieved profile is a weighted sum of the a priori and the true profile, with the weights provided by the AVK. When the information content of the measurement is low or the retrieval is too heavily constrained, the AVK will tend toward zero and the solution will remain close to the a priori. Conversely, when the information content is high or the retrieval loosely constrained, the AVK will approach the unit matrix. The optimal estimation AVK is therefore a measure of how much information is extracted from the measurement, with its trace commonly denoted “degrees of freedom for signal”. By contrast, the total column averaging kernels introduced above are not a measure of how much information is extracted from the measurement, as they accompany an unconstrained retrieval. A perfect measurement would correspond to an all-ones vector A_z=1 for all z. However, inherent to (infrared) sounding, sensitivity varies as a function of thermal contrast, and thus altitude, so this ideal can never be met.

The second important difference relates to the fact that a vertical profile is retrieved in the optimal estimation type retrievals. The rows of A express the vertical resolution of the retrieval, where the ideal is a narrow function that peaks at its corresponding altitude. As no vertical information is extracted in the total column retrievals, this interpretation also does not apply. A_z is a (normalized) vector that shows which layers of the atmosphere offer in principle the greatest sensitivity. It naturally peaks high up in the atmosphere, irrespective of the location of the trace gas. However, combining the total column AVK with the a priori profile does allow extracting a vector which characterizes the vertical sensitivity, as we now show. Starting from Eqs. (26) and (28) we have

$$\begin{array}{}\text{(34)}& {\displaystyle }\text{HRI}={\text{SF}}^a\left({\widehat{X}}^a-B\right)\text{(35)}& {\displaystyle }=\sum _z{\text{SF}}_z\left({\widehat{X}}^a-B\right)a_z\text{(36)}& {\displaystyle }=\sum _z{\text{SF}}_z\left({\widehat{X}}_z^a-B_z\right)\text{(37)}& {\displaystyle }=\sum _z{\widehat{\text{HRI}}}_z^a,\end{array}$$

where we defined ${\widehat{\text{HRI}}}_z^a={\text{SF}}_z\left({\widehat{X}}_z^a-B_z\right)$ similarly to Eq. (24) where we defined ${\text{HRI}}_z={\text{SF}}_z(X_z-B_z)$. A retrieval that assumes a given a priori vertical profile a_z therefore implicitly assumes that the HRI (the trace gas signal) can be decomposed into partial ${\widehat{\text{HRI}}}_z^a$ corresponding to spectral change at each altitude. Note that while $\text{HRI}={\sum }_z{\widehat{\text{HRI}}}_z^a={\sum }_z{\text{HRI}}_z$, the same does not necessarily hold for each individual level ${\widehat{\text{HRI}}}_z^a\ne {\text{HRI}}_z$, as the assumed profile can differ from the actual profile. With this we can define the normalized assumed HRI profile or, equivalently, the probable vertical partitioning $V_z^a$ of the signal:

$$\begin{array}{}\text{(38)}& V_z^a={\displaystyle \frac{{\widehat{\text{HRI}}}_z^a}{\text{HRI}}},={\displaystyle \frac{{\text{SF}}_z\left({\widehat{X}}_z^a-B_z\right)}{{\text{SF}}^a({\widehat{X}}^a-B)}}=A_z^a{a}_z.\end{array}$$

From the first equality, it follows that ${\sum }_z{V}_z^a=\mathrm{1}$ (and also follows the last term and Eq. 29). The three profiles a_z, $A_z^a$ and $V_z^a$ are illustrated in Fig. 2 for a typical NH₃ retrieval. As can be seen, the maximum of $V_z^a$ is shifted upwards compared with a_z, due to the more favourable thermal contrast higher up in the atmosphere.

3.4 AVK application

There are two alternative ways in which averaging kernels can be exploited to remove the impact of the vertical profile assumption of the retrieval (Palmer et al., 2001; Eskes and Boersma, 2003; Cooper et al., 2020).

3.4.1 Method 1: simulating measurements of the modelled columns

Let $M_z^m$ be a modelled profile with corresponding total column M^m and normalized profile (enhancement):

$$\begin{array}{}\text{(39)}& m_z={\displaystyle \frac{M_z^m-B_z}{M^m-B}}.\end{array}$$

We can simulate what would have been retrieved if the modelled profiles were observed using Eq. (32):

$$\begin{array}{}\text{(40)}& {\widehat{M}}^a=\sum _z{A}_z^a(M_z^m-B_z)+B.\end{array}$$

This ${\widehat{M}}^a$ is directly comparable with ${\widehat{X}}^a$ as the same a priori profile shape a_z is used for both retrievals. However, in case the a priori profile significantly differs from the truth, both ${\widehat{X}}^a$ and ${\widehat{M}}^a$ can deviate far from the truth.

3.4.2 Method 2: using modelled vertical profiles as a priori

Rather than altering the modelled column, an attractive alternative is to alter the retrieved column to use instead of the a priori vertical profile, the modelled profile (see Boersma et al., 2016, Appendix D):

$$\begin{array}{}\text{(41)}& {\displaystyle }{\widehat{X}}^m={\displaystyle \frac{\text{HRI}}{{\text{SF}}^m}}+B\text{(42)}& {\displaystyle }={\displaystyle \frac{{\text{SF}}^a({\widehat{X}}^a-B)}{{\sum }_z{\text{SF}}_z{m}_z}}+B\text{(43)}& {\displaystyle }={\displaystyle \frac{{\widehat{X}}^a-B}{{\sum }_z{A}_z^a{m}_z}}+B.\end{array}$$

The averaging kernel associated with ${\widehat{X}}^m$ is

$$\begin{array}{}\text{(44)}& A_z^m={\displaystyle \frac{{\text{SF}}_z}{{\text{SF}}^m}}=A_z^a{\displaystyle \frac{{\text{SF}}^a}{{\text{SF}}^m}}.\end{array}$$

This ${\widehat{X}}^m$ can be directly compared with M^m, as both employ the same profile m_z. Note that ${\widehat{X}}^m$ depends only on the shape of the modelled profile, not the total column. This method can be used to obtain an improved retrieval by using a modelled profile that approaches the reality better than a static a priori profile. Since it was first applied to NH₃ (Whitburn et al., 2016), the ANNI retrieval has been capable of using modelled profiles by adapting the input parameters to the network. However, when the modelled profiles were changed, the entire retrieval process had to be redone. Using Eq. (43) and the provided AVKs, changing the a priori profile can be done a posteriori by the data users. An important practical note is that the post-filter of the retrieval (see Sect. 2) includes a threshold on the scaling factor, and that this filter should be reevaluated for ${\widehat{X}}^m$ using SF^m.

Both methods can be summarized as

$$\begin{array}{}\text{(45)}& {\displaystyle }{M}^m\stackrel{\text{method}\mathrm{1}}{⟶}{\widehat{M}}^a\text{to be compared with}{\widehat{X}}^a\text{(46)}& {\displaystyle }{\widehat{X}}^a\stackrel{\text{method}\mathrm{2}}{⟶}{\widehat{X}}^m\text{to be compared with}{M}^m.\end{array}$$

Equations (39), (40) and (43) can be combined as

$$\begin{array}{}\text{(47)}& {\displaystyle }{\widehat{X}}^m={\displaystyle \frac{({\widehat{X}}^a-B)(M^m-B)}{{\sum }_z{A}_z^a(M_z^m-B_z)}}+B\text{(48)}& {\displaystyle }={\displaystyle \frac{({\widehat{X}}^a-B)(M^m-B)}{{\widehat{M}}^a-B}}+B,\end{array}$$

or

$$\begin{array}{}\text{(49)}& {\displaystyle \frac{{\widehat{X}}^m-B}{M^m-B}}={\displaystyle \frac{{\widehat{X}}^a-B}{{\widehat{M}}^a-B}},\end{array}$$

which shows that, when the goal is to compare the ratio between model and retrieved columns, both methods are equivalent (Cooper et al., 2020).

Figure 3 AVK normalization factors (N) for 1 day (17 June 2015, a) and for 1 year (2015, b) of IASI observations. In both cases, data originates from the morning overpass of IASI-A. The insets show the respective normalized histograms.

3.5 Practical considerations

In the ANNI retrieval formalism, the total scaling factor SF^a is calculated directly by the NN, not via intermediate SF_z and application of Eq. (28). However, these SF_z, which are required to calculate the AVKs, can be estimated by exploiting the flexibility of the NN. For NH₃, the NN is trained for a wide variety of Gaussian profiles, with peak altitudes ranging from 0 to 20 km and σ from 100 m to 3 km. The SF_z can be estimated from the network using the input parameters z_peak=z and σ=100 m for the Gaussian profile. For this calculation, an HRI input parameter is also required and the choice was made to use the observed HRI. The corresponding column that can be calculated from this satisfies

$$\begin{array}{}\text{(50)}& {\widehat{X}}^{\mathrm{|}z}={\displaystyle \frac{\text{HRI}}{{\text{SF}}^{\mathrm{|}z}}}+B,\end{array}$$

where ${\widehat{X}}^{\mathrm{|}z}$ is the retrieved total column assuming all the trace gas enhancement is situated in the narrow Gaussian band around altitude z. SF^|z is the corresponding total scaling factor, which is used to approximate the local scaling factor SF_z. With this the AVK can be constructed as

$$\begin{array}{}\text{(51)}& A_z^a={\displaystyle \frac{{\text{SF}}_z}{{\text{SF}}^a}}\approx {\displaystyle \frac{{\text{SF}}^{\mathrm{|}z}}{{\text{SF}}^a}}={\displaystyle \frac{({\widehat{X}}^a-B)\text{HRI}}{({\widehat{X}}^{\mathrm{|}z}-B)\text{HRI}}}={\displaystyle \frac{{\widehat{X}}^a-B}{{\widehat{X}}^{\mathrm{|}z}-B}}.\end{array}$$

The formulas provided above are exact in the linear limit, but for large columns, SF^a and SF^|z have an increasingly high dependence on the value of the HRI. The NN takes into account this dependence so that Eqs. (2) and (50) are always good approximations of the true X, provided that either the assumed a priori profile is correct or that the tracer is confined to a narrow layer. However, there is no guarantee that

$$\begin{array}{}\text{(52)}& {\text{SF}}^a\stackrel{\mathrm{?}}{=}\sum _z{\text{SF}}_z{a}_z,\end{array}$$

or thus that ${\sum }_z{A}_z^a{a}_z$ equals one. A consistent AVK can be obtained, however, as

$$\begin{array}{}\text{(53)}& A_z^a={\displaystyle \frac{\mathrm{1}}{N}}{\displaystyle \frac{{\widehat{X}}^a-B}{{\widehat{X}}^{\mathrm{|}z}-B}},\end{array}$$

with $N={\sum }_z{A}_z^a{a}_z$ a normalization factor. This factor also guarantees that applying the averaging kernels on the a priori vertical profile returns the retrieved column:

$$\begin{array}{}\text{(54)}& \sum _z{A}_z^a\left({\widehat{X}}_z^a-B_z\right)+B=\sum _z{A}_z^a{a}_z\left({\widehat{X}}^a-B\right)+B={\widehat{X}}^a.\end{array}$$

The normalization factors are shown in Fig. 3 for 1 day and 1 year of global measurements. The histograms shown in the insets indicate average values around 0.98–0.99 (±0.03–0.04), illustrating the consistency of the approach and the fact that non-linear effects are modest. The areas where the normalization factors are furthest from 1 are affected by low clouds (e.g. off the west coast of South Africa) or sea ice.

The necessity of using normalization factors follows from the fact that in the non-linear regime, the SF_z are not uniquely defined, depending on the concentration in the layer z (non-linearity) and the other layers (non-additivity). However, as we have shown, in the neighbourhood of the solution, a fully consistent AVK can be obtained after renormalization. It is this AVK that is recommended when applying method 1 in model comparisons. However, in case method 2 is employed and when the modelled profile concerns a narrow layer at high altitude (e.g. a pyro-convective fire plume) or significantly deviates from the a priori, it can be better not to renormalize. In particular, when we have a narrow modelled profile layer at an altitude z^′, with $m_z={\mathit{\delta }}_{z{z}^{\prime }}$, the second method, without renormalization, yields the expected

$$\begin{array}{}\text{(55)}& {\displaystyle }{\widehat{X}}^m={\displaystyle \frac{{\widehat{X}}^a-B}{{\sum }_z{A}_z^a{\mathit{\delta }}_{z{z}^{\prime }}}}+B\text{(56)}& {\displaystyle }={\displaystyle \frac{{\widehat{X}}^a-B}{A_{z^{\prime }}^a}}+B\text{(57)}& {\displaystyle }=\left({\widehat{X}}^{\mathrm{|}{z}^{\prime }}-B\right)+B\text{(58)}& {\displaystyle }={\widehat{X}}^{\mathrm{|}{z}^{\prime }}.\end{array}$$

The output files of the ANNI retrieval contain ${\widehat{X}}^a$, ${\widehat{X}}^{\mathrm{|}z}$, B_z and N. With this, $A_z^a$ can be calculated if required from Eq. (51) or renormalized via Eq. (53).

4 Temporal consistency

4.1 Pseudoinverse

The generalized error covariance matrix S plays a key role in the calculation of the HRI. As a symmetric matrix, S has real eigenvalues λ_i and can be decomposed as

$$\begin{array}{}\text{(59)}& \mathbf{S}=\sum _{i=\mathrm{1}}^n{\mathit{\lambda }}_i{\mathit{s}}_{\mathit{i}}{{\mathit{s}}_{\mathit{i}}}^{\mathrm{T}},\end{array}$$

with all s_i orthogonal to each other. It follows that its inverse can be written as

$$\begin{array}{}\text{(60)}& {\mathbf{S}}^{-\mathrm{1}}=\sum _{i=\mathrm{1}}^n{\displaystyle \frac{\mathrm{1}}{{\mathit{\lambda }}_i}}{\mathit{s}}_{\mathit{i}}{{\mathit{s}}_{\mathit{i}}}^{\mathrm{T}}.\end{array}$$

This expression leads to an intuitive interpretation of the HRI (Clarisse et al., 2013): it can be seen as a weighted projection of the spectrum onto the Jacobian, with the directions that usually exhibit the most variability carrying the lowest weight.

The distribution of the eigenvalues of the covariance matrix used for the NH₃ HRI (with a spectral range covering 812–1126 cm⁻¹) is shown in Fig. 4. Three domains can be distinguished: (i) the 30 highest values corresponding to the principal components, (ii) around 1200 values corresponding mainly to instrumental noise and (iii) around 20 very small eigenvalues. These smallest eigenvalues are of the order of the numerical precision at which the covariance matrix is calculated and, in essence, correspond to directions not occurring in IASI spectra. While random instrumental noise would be expected to occur in all directions, apodization and L1 post-processing remove some. Such directions carry the most weight in S⁻¹, but as they are not found in real spectra they do not contribute much to the total HRI (as can easily be verified numerically).

Figure 4 Eigenvalue spectrum of the covariance matrix S used for the calculation of the HRI of NH₃. The eigenvalues are ordered from largest to smallest.

However, small changes to the instrument calibration or post-processing can alter the contribution of these directions in the IASI spectra, and because they carry such a large weight in the HRI, they can affect its value considerably. This explains why the HRI in the past has been found to be very sensitive to such changes (Van Damme et al., 2021). It also explains the occurrence of (small) biases between the different instruments. The solution is fortunately simple (Rodgers and Connor, 2003; Eaton, 2007) and is obtained by disregarding the terms corresponding to the very small eigenvalues in Eq. (60). As we will show later, using such a pseudoinverse does not eliminate the effects of L1C changes completely but reduces their magnitude considerably.

After the pseudoinverse was implemented, an unexpected change was observed in the value of the HRIs on spectra from the period on which the covariance matrix was built. It turns out that while the scalar product of observed IASI spectra with the eigenvectors corresponding to the lowest eigenvalues ${\mathit{L}}^T{\mathit{s}}_{\mathit{i}}/{\mathit{\lambda }}_i$ is near zero, this is not the case for spectra generated with the forward model. This is due to small discrepancies between spectra generated by the forward and actual spectra that are magnified by the $\mathrm{1}/{\mathit{\lambda }}_i$ factors. Hence, synthetic HRIs calculated on the training set of the neural network have in the past been overestimated, resulting in low biases in the retrieved columns. As we will show in Sect. 6, the magnitude of this bias was around 18 % in ANNI-v3. For the retrieval of other trace gases presented in Franco et al. (2018), especially those operating on a smaller spectral range, the bias has been evaluated to be much smaller.

4.2 Carbon dioxide

As the mean spectrum and covariance matrix that are used for the HRI are calculated from spectra measured within 1 reference year (2013), long-term changes in atmospheric composition that affect the spectral region of interest can have unwanted effects on the HRI. This was first noted in Van Damme et al. (2021), where a spurious trend was seen in the HRI NH₃ data over remote regions. It was attributed to the increase in global carbon dioxide (CO₂) concentrations, because of the presence of a weak CO₂ absorption band in the 920–990 cm⁻¹ spectral region (Whitburn et al., 2021) where NH₃ has its strongest absorption. A linear correction on the HRI of the order of 0.03 yr⁻¹ was introduced to compensate for this effect. However, because of seasonal variations, and possible temperature dependence of the interference an HRI which is less sensitive to CO₂ changes is preferable. One option is to build the covariance matrix from spectra spanning the entire period of IASI measurements.

An alternative approach is to account directly for the effects of CO₂ in the calculation of the HRI. The HRI formula is related to generalized least squares estimation and can be expanded to include multiple variables that are simultaneously estimated (Walker et al., 2011; Theys et al., 2022). In our case, the Jacobian vector becomes a two-column matrix, one column corresponding to NH₃ and the other to CO₂. The HRI formula Eq. (1) remains formally identical (with only the first component of the two-element HRI vector of interest). The effect of this change on the long-term trend of the HRI is detailed in Sect. 4.4.

Figure 5 Monthly average HRI time series over 10 remote regions for the three IASI instruments separately. The top panel shows the uncorrected time series, and the other panels, from top to bottom, show the effects of the corrections that are applied consecutively.

4.3 Cloud clearing

The ERA5 model output replaces satisfactorily the IASI L2 for all input parameters, except for the surface temperature and cloud cover. These are spatially and temporally too variable for model output to be representative for an IASI footprint at a given time. All previous reanalysed ANNI-NH₃ products still relied on different versions of the IASI L2 cloud product. Recently, Whitburn et al. (2022) developed an NN-based cloud flag. Trained with data from the latest version (v6.5) of the official L2 cloud product, it inherits all its advantages as a proven and well-validated product. The NN utilizes carefully selected IASI channels as input (excluding channels affected by long-lived tracers CO₂, N₂O, CH₄, CFC-11 and CFC-12) and was shown to be temporally consistent, and coherent across the three IASI instruments. The network presented in Whitburn et al. (2022) was trained to distinguish completely clear scenes (0 % cloud cover) from the rest. For the NH₃ processing, two additional networks were trained to distinguish scenes with a cloud cover below 10 % and 25 % respectively. With this, three cloud flags are available and these have now been integrated in v4 of the reanalysed ANNI-NH₃ product. The results presented in the rest of this paper utilize the flag corresponding to the 10 % threshold.

4.4 Residual bias corrections

The stability of the HRI was evaluated over 10 remote regions where only background columns of NH₃ are expected. Their average monthly HRI is shown in Fig. 5 for the three IASI instruments separately. The top panel shows the average as obtained with the HRI setup as described above, i.e. with pseudoinverse and with a CO₂ Jacobian. As with previous versions of the product, a spurious linear trend is observed, but thanks to the introduction of the CO₂ Jacobian, its magnitude is reduced to about 0.01 yr⁻¹, compared with 0.03 yr⁻¹ previously. A slightly steeper decrease is observed for Metop B. We correct for these trends by adding a time-dependent offset as in Van Damme et al. (2021). The result after correction is shown in the second panel of Fig. 5.

A detailed analysis was made of this time series to detect offsets between the different instruments and shifts that coincide with known changes in the IASI L1C data. The largest of these shifts is the offset of 0.11 seen between IASI-C and the two other instruments. Small offsets in the HRI time series of IASI-A were found in 2010, 2015 and 2017, and in the HRI of IASI-B in 2015. For each of these, offset corrections were calculated in the range of 0.01–0.03. Thanks to the pseudoinverse, their magnitude is drastically reduced. Previously, offsets as large as 0.6 were observed. The resulting corrected time series is shown in the third panel. This time series is temporally stable and shows an excellent consistency between the three instruments, but it exhibits a weak seasonal cycle, likely due to the combined effect of seasonal changes in the concentrations of H₂O and volatile organic compounds that absorb in the same spectral range as NH₃. To remove this seasonality an offset depending on latitude and month of the year was calculated from 2012 to 2014 IASI-A data and applied on all data. The HRI after correction is shown in the bottom panel of Fig. 5. Thanks to the improved setup of the HRI, and the new cloud product, the magnitude and therefore also uncertainty of all these corrections is lower than in the previous product, which results in much improved temporal consistency.

Figure 6 An NH₃ seasonal average derived from $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ monthly averages of the reanalysis product of ANNI v4. Data include all measurements from IASI-A (October 2007 to December 2019), IASI-B (March 2013 to September 2022) and IASI-C (September 2019 to September 2022), with a cloud fraction below 10 %.

5 Other changes to the retrieval network

An additional change in the setup of the HRI concerns the choice of the spectra used for determining the mean background spectra and its associated covariance matrix. As before (Franco et al., 2018) we use a random selection of IASI spectra from the year 2013, but now with a proportionally larger number of spectra from selected parts of the Sahara, Arabian, Great Australian and Namib deserts. It was found that this was an efficient way for countering the small negative biases that are seen over these areas and associated with surface emissivity variations. It also leads to a better detection of NH₃ transport over deserts.

Since ANNI v2 (Van Damme et al., 2017), the reanalysis product relies on a surface temperature retrieved from a custom-built neural network, rather than the IASI L2 surface temperature. With ANNI v4, this network has been retrained from data that were generated using the latest version (v6.5) of the IASI L2 algorithm. The input parameters of the NN for the retrieval of surface temperature include 60 selected baseline channels (a subset of the channels used in the cloud NN), surface altitude, total water vapour column and the three output values of the cloud NNs. Mean and standard deviation of the difference between the L2 surface temperature and that retrieved from the network are of the order of 0.5 and 1.5 K respectively for cloud fraction up to 25 %.

A final series of changes concern the network architecture and training database. In previous versions, separate neural networks were employed for the retrieval over land and ocean. These networks were trained respectively with Gaussian a priori profiles peaking at the surface, and a more general one, with profiles peaking at different altitudes. However, a careful comparison showed that the more general network performed as well as the network trained specifically for profiles peaking at the surface. For this reason, only one network was trained for version 4, for a priori profiles peaking at altitudes z₀ from 0 to 20 km, with a width σ in the range of 0.1–3 km. In view of the averaging kernel calculation, 20 % of the profiles of the training database have an NH₃ profile with a σ of 100 m.

6 Evaluation

6.1 Comparison with version 3

As outlined before, v4 has an improved temporal consistency compared with v3. In this section, we provide a short assessment of the new NH₃ spatial distributions and how they compare with previous versions. As an illustration of the new product, a seasonal average over 2007–2022 is presented in Fig. 6. The distributions follow closely the ones of previous versions (Van Damme et al., 2015). Comparisons with version 3 are provided in Figs. 7 and 8. The main differences are: (1) overall larger columns, especially in areas with high columns. As explained in Sect. 4.1, this is due to the more robust way of calculating the HRI, which makes it less sensitive to small errors in the forward model. Comparing individual observations, the new version is about 10 % (for low HRI) to 20 % (for high HRI) larger. (2) A notable improvement over high latitudes thanks to the improved bias correction (see Sect. 4.4). Averaged columns were clearly overestimated in ANNI v3 for such observations, especially over Greenland and Antarctica but also Canada and Russia. (3) Slightly higher concentrations over deserts, partly due to the overall increase in v4 and partly due to the larger weight of deserts in the construction of mean background spectra and associated covariance matrix (see Sect. 5). In v3, negative HRIs were consistently observed over certain deserts, resulting in negative average columns or low biases. This problem is not entirely gone (e.g. the average HRI over parts of the Arabian desert is still negative), but much improved, resulting in more consistent pronounced transport patterns over, for example, the Sahara Desert. (4) Overall larger concentrations over oceans. The improved bias correction (last step of the HRI correction presented in Sect. 4.4) enables removal of practically all negative values on a long-term average. Average columns over remote ocean are now around 7×10¹⁴ molec. cm⁻² almost everywhere. An adverse effect of the correction might be overestimated columns over the Red Sea, Persian Gulf and Mediterranean Sea.

Figure 7 Comparison between retrieved columns with ANNI v3 and v4 for all morning observations of 17 June 2015.

Figure 8 Comparison between the NH₃ columns of the near-real-time products of ANNI v3 (a) and v4 (b) on a $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ grid. The inset shows the difference between the two maps. Data include all morning IASI-A data from 2008–2018, with a cloud fraction below 10 %. Parallels are drawn every 15^∘ and meridians every 30^∘.

The most obvious remaining artefact in the v4 distribution concerns the continuity of the land–sea transitions. While they are reasonable for some regions of outflow (Gulf of Mexico, Mediterranean Sea), off the west coast of Africa, over the Arabian Sea, Gulf of Bengal or Yellow Sea, the transition is too abrupt to be realistic. The origin of this problem is that different NH₃ profiles are used for land and ocean. With the introduction of the AVKs, this does not constitute a problem in model comparison or assimilation. However, for stand-alone use of the product, it would be desirable in the future to improve the parametrization of the a priori vertical profile shape.

6.2 Comparison with an optimal estimation retrieval

Given the low bias in ANNI NH₃ v3, it is important to exclude the presence of other biases related to the HRI calculation. Here, we present the results of an independent intercomparison that was conducted between the ANNI v4 retrieval output and that of an optimal estimation approach which relies on spectral fitting.

For the optimal estimation retrieval, the Atmosphit forward and inverse model was used (Coheur et al., 2005), which is the same tool whose forward model is used for the construction of the ANNI training database. The optimal estimation was set up as follows: the retrieval range was set to 900–975 cm⁻¹. Total columns of NH₃ were retrieved with a fixed vertical profile, using the same parametrization as in ANNI NRT. The NH₃ variance was set to 1000 %, corresponding to an almost unconstrained retrieval. Together with NH₃, H₂O was retrieved in 10 partial columns, with the a priori coming from the IASI L2. Total columns of CO₂, O₃ and CFC-12 were retrieved as well as the surface temperature. Spectral emissivity was taken from Zhou et al. (2013). Before presenting the results, it should be emphasized that despite the similarities in both retrieval approaches (same input parameters, vertical profiles, forward model), no perfect agreement is expected because of (1) use of a narrower spectral range in the optimal estimation retrieval, (2) different propagation of instrumental noise to the retrieval result, (3) limitations of the fitting model (e.g. with respect to fitting water vapour or surface emissivity) and (4) errors related to the imperfect training of the neural network.

For the comparison, 2 days were selected, one over Europe and one over North America, with relatively high NH₃ columns. The results are shown in Fig. 9. Intercepts, mean and median differences are all of the order of 10¹⁵ molec. cm⁻² or below. Regression slopes, calculated using iteratively reweighted least squares to remove the impact of outliers, are 0.97 and 1.05. While the scatter around the 1–1 lines is not negligible (with standard deviation of the differences around 3–4 ×10¹⁵ molec. cm⁻²), these numbers demonstrate the overall consistency of both retrieval approaches and do not indicate a significant bias. Note that the ANNI retrieval approach has numerous advantages over optimal estimation, as discussed in Whitburn et al. (2016).

The last detailed global validation of the IASI NH₃ product was based on a comparison of ground-based FTIR measurements of NH₃ with the LUT-based NH₃ product, where a low bias around 35 % was found (Dammers et al., 2016). Since then, two independent validation studies have been conducted. One study (Guo et al., 2021) compared IASI ANNI v3 with in situ measurements in Colorado, USA, and found regression slopes ranging from 0.78 to 1.1 and intercepts of the order of 1–2×10¹⁵ molec. cm⁻². A second study (Wang et al., 2020) compared IASI NH₃ columns with columns obtained from FTIR measurements in Hefei, China. Here, mean differences around 3.5×10¹⁵ molec. cm⁻² (IASI being lower) were found and regression slopes close to one. Given the results of the comparison with the optimal estimation method, we do not expect any significant bias in v4 for columns above 1×10¹⁶ molec. cm⁻² in comparisons that correct for the vertical profile assumption of the retrieval. A comprehensive validation of the v4 product is foreseen within the framework of ESA's CCI+ precursors for the aerosol and ozone ECV project, which should confirm this, as well as assess the performance of the algorithm on low columns. Apart from validation of the columns in an absolute or relative sense, comparison with FTIR columns will also allow evaluating regional NH₃ trends derived from IASI data. Such an evaluation could also be made with bottom-up inventories or with data derived from in situ measured concentrations. However, in that case, there is the additional difficulty that long-term trends of other inorganic pollutants (NO₂, SO₂) affect NH₃ columns differently than emissions or local concentrations (e.g. Lachatre et al., 2019), necessitating the intervention of a (chemistry transport) model.

7 Uncertainties

7.1 Propagation of uncertainty

In previous ANNI versions, an estimated uncertainty ${\mathit{\sigma }}_{\widehat{X}}$ was calculated for each individual measurement $\widehat{X}$ via (Ku, 1966)

$$\begin{array}{}\text{(61)}& {\mathit{\sigma }}_{\widehat{X}}^{\mathrm{2}}=\sum _i{\left({\displaystyle \frac{\partial \widehat{X}}{\partial p_i}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_{p_i}^{\mathrm{2}},\end{array}$$

with ${\mathit{\sigma }}_{p_i}$ the uncertainties of the different input parameters p_i. This equation assumes uncorrelated uncertainties, but as this cannot always be justified, in ANNI v4, we switch to the more general (Tellinghuisen, 2001)

$$\begin{array}{}\text{(62)}& {\mathit{\sigma }}_{\widehat{X}}^{\mathrm{2}}={\mathit{J}}^{\mathrm{T}}{\mathbf{S}}_p\mathit{J},\end{array}$$

with S_p the error covariance matrix of the input parameters (with covariances $S_{p,ij}={\mathit{\sigma }}_{p_i{p}_j}$) and J the Jacobian of the retrieval, with components $\frac{\partial \widehat{X}}{\partial p_i}$.

In the ANNI retrieval framework, the input parameters include the skin temperature, the surface pressure, the HRI, the surface emissivity, the zenith angle, the width and the peak of the Gaussian vertical NH₃ profile, the temperature profile (15 levels) and the water vapour profile (7 levels). After some preliminary analysis, it was concluded that only the correlations between the uncertainties in the temperature profile cannot be neglected. We therefore employ a block diagonal covariance matrix, block diagonal for the elements pertaining to the temperature profile, and diagonal for all other input parameters. As for uncertainty on the vertical profile, this source of uncertainty is removed when applying averaging kernels. For this reason, uncertainties are reported with and without the vertical profile uncertainty, to be used according to whether AVKs are applied.

Figure 9 Comparison between ANNI v4 NH₃ columns and retrievals based on optimal estimation for two scenes, one over Europe (a–c: 18 April 2013, Metop A morning overpass) and North America (d–f: 6 May 2021, Metop B morning overpass). Panels (a) and (d) depict the optimal estimation retrieved columns. Panels (b) and (e) are scatter plots between the two retrievals, where each observation is colour coded according to thermal contrast (brightness temperature of the surface minus the temperature at half the boundary layer height). Panels (c) and (f) summarize the comparison by means of histograms of the differences.

7.2 Random and systematic uncertainties

In total, we report four types of uncertainty for each observation: random or systematic, and with or without the vertical profile uncertainty included. Reporting random and systematic uncertainties separately is a generally recommended practice (Boersma et al., 2004; Merchant et al., 2017; Sayer et al., 2020). Random uncertainties describe errors specific to a single measurement, and assuming a normal distribution, these average out over many repeated measurements. Systematic uncertainties are those that exhibit correlations in time or space, and are thus associated with more than one measurement. This type of error can lead to biases in the measurement dataset. In ANNI v4, both random ${\mathit{\sigma }}_{\mathrm{r}\widehat{X}}$ and systematic ${\mathit{\sigma }}_{\mathrm{s}\widehat{X}}$ uncertainties are calculated using Eq. (62) and estimates of the random and systematic uncertainties and covariances of the input parameters.

Random and systematic uncertainties can be combined and averaged in different ways according to the needs of the user. In particular, for a given measurement $\widehat{X}$, a total uncertainty estimate can be obtained as (Gomez-Pelaez et al., 2013)

$$\begin{array}{}\text{(63)}& {\mathit{\sigma }}_{\widehat{X}}^{\mathrm{2}}={\mathit{\sigma }}_{\mathrm{r}\widehat{X}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{s}\widehat{X}}^{\mathrm{2}}.\end{array}$$

An average measurement uncertainty can be associated with an average $\stackrel{\mathrm{‾}}{X}$ of a series of n measurements ${\widehat{X}}_i$ as

$$\begin{array}{}\text{(64)}& {\displaystyle }{\mathit{\sigma }}_{\stackrel{\mathrm{‾}}{X}}^{\mathrm{2}}={\mathit{\sigma }}_{\mathrm{r}\stackrel{\mathrm{‾}}{X}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{s}\stackrel{\mathrm{‾}}{X}}^{\mathrm{2}}\text{(65)}& {\displaystyle }=\sum _{i=\mathrm{1}}^n{\left({\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{r}{\widehat{X}}_i}}{n}}\right)}^{\mathrm{2}}+{\left(\sum _{i=\mathrm{1}}^n{\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{s}{\widehat{X}}_i}}{n}}\right)}^{\mathrm{2}}.\end{array}$$

For the special case where all random and systematic uncertainties are the same, we obtain

$$\begin{array}{}\text{(66)}& {\mathit{\sigma }}_{\stackrel{\mathrm{‾}}{X}}^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{n}}{\mathit{\sigma }}_{\mathrm{r}\stackrel{\mathrm{‾}}{X}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{s}\stackrel{\mathrm{‾}}{X}}^{\mathrm{2}},\end{array}$$

which tends to the expected ${\mathit{\sigma }}_{\stackrel{\mathrm{‾}}{X}}={\mathit{\sigma }}_{\mathrm{s}\stackrel{\mathrm{‾}}{X}}$ for large n.

7.3 Uncertainties of the input parameters

As most input parameters come without an uncertainty budget, let alone covariances, we made best-effort estimates of the co(variance) based on the limited information that is available. For now, the same (co)variances were used for the near-real time as for the reanalysed NH₃ product. It is also important to note that the systematic uncertainties of the input parameters vary according to the time and space scales that are considered (Boersma et al., 2004; Merchant et al., 2017; Sayer et al., 2020). Temperature profiles, for example, may be more biased monthly than annually. Here, we estimate systematic uncertainties with a typical L3 gridded data product in mind, i.e. for spatial scales of the order of 1^∘ latitude and longitude or less, and for time periods of the order of 1 month or less.

The (co)variances, summarized in Table 2, were determined as follows:

HRI

By definition, the random uncertainty on the HRI equals one. We estimate a systematic uncertainty of 0.1 due to potential and residual interferences (e.g. surface emissivity, other species). To account for potential biases in the spectroscopy and generalized error covariance matrix, we add to this an additional 10 % on the calculated HRI value.

Skin temperature

Random and systematic uncertainties were set to 1.5 and 0.5 K respectively. These values are in line with the difference between the IASI L2 skin temperature product and the dedicated neural network used for the reanalysis product of ANNI (see Sect. 5).

Emissivity

For emissivity, which originates from the monthly climatology of Zhou et al. (2013), an uncertainty of 0.01 and 0.005 was assumed for respectively the random and systematic components.

Temperature profile

Variances were set based on validation results of the IASI level 2 (EUMETSAT, 2021): systematic uncertainties of 1 K for the surface level and 0.5 K for the other levels; random uncertainties of 2 K for the surface level and 1 K for the other levels for land observations, and 1 K for the surface level and 0.5 K for the other levels for ocean observations. Covariance matrices were then built by appropriate scaling of correlation matrices. These were built based on a statistical analysis of the differences between collocated ERA5 and IASI L2 profiles. Correlation coefficients were set to 0.5 between neighbouring levels and 0.25 between levels that are two levels apart. Above 10 km, no strong correlations were observed, and the covariance was therefore assumed to be diagonal for these levels.

Water vapour profiles

Relying again on the IASI level 2 validation report (EUMETSAT, 2021), random uncertainties were set to 10 % below 3 km and 20 % above. Systematic uncertainties were set to half these numbers.

Surface pressure

A random and systematic uncertainty of 500 and 250 Pa was used.

NH₃ profiles

The uncertainties related to the NH₃ profile arise from uncertainties associated with the width and the peak of the Gaussian shaped vertical profile. Random and systematic uncertainties of 200 and 100 m were used for both parameters. Given the short lifetime of NH₃ in the atmosphere, these are likely of the right order of magnitude. To obtain better estimates in the future, a thorough analysis using in situ measurements or modelled profiles would be desirable.

Table 2 Estimated random and systematic uncertainties of the input parameters.

7.4 Uncertainty budget of NH₃

It is useful, remembering the general form ${\widehat{X}}^a=\text{HRI}/{\text{SF}}^a+B$ of the retrieval, to rewrite the propagation of uncertainty in terms of the uncertainty of the nominator and denominator (see also Boersma et al., 2004; van Geffen et al., 2022). Neglecting the small dependence of the SF on the HRI, we obtain

$$\begin{array}{}\text{(67)}& {\displaystyle }{\mathit{\sigma }}_{\widehat{X}}^{\mathrm{2}}={\left({\displaystyle \frac{\partial \widehat{X}}{\partial \text{HRI}}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_{\text{HRI}}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial \widehat{X}}{\partial \text{SF}}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_{\text{SF}}^{\mathrm{2}}\text{(68)}& {\displaystyle }={\displaystyle \frac{{\mathit{\sigma }}_{\text{HRI}}^{\mathrm{2}}}{{\text{SF}}^{\mathrm{2}}}}+{\displaystyle \frac{{\mathit{\sigma }}_{\text{SF}}^{\mathrm{2}}}{{\text{SF}}^{\mathrm{2}}}}(X-B{)}^{\mathrm{2}}.\end{array}$$

Taking into account both random and systematic uncertainties, we see from Table 2 that the uncertainty on the HRI has an absolute (constant) and a relative (proportional to the value of the HRI) component, so that

$$\begin{array}{}\text{(69)}& {\displaystyle }{\mathit{\sigma }}_{\widehat{X}}^{\mathrm{2}}={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{abs},\mathrm{HRI}}^{\mathrm{2}}}{{\text{SF}}^{\mathrm{2}}}}+{\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{rel},\mathrm{HRI}}^{\mathrm{2}}}{{\text{SF}}^{\mathrm{2}}}}+{\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{SF}}^{\mathrm{2}}}{{\text{SF}}^{\mathrm{2}}}}(X-B{)}^{\mathrm{2}}\text{(70)}& {\displaystyle }=\underset{{\mathit{\sigma }}_{\mathrm{abs}}^{\mathrm{2}}}{\underbrace{{\displaystyle \frac{{\mathrm{1}}^{\mathrm{2}}+{\mathrm{0.1}}^{\mathrm{2}}}{{\text{SF}}^{\mathrm{2}}}}}}+\underset{{\mathit{\sigma }}_{\text{rel}}^{\mathrm{2}}}{\underbrace{\left({\mathrm{0.1}}^{\mathrm{2}}+{\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{SF}}^{\mathrm{2}}}{{\text{SF}}^{\mathrm{2}}}}\right)(X-B{)}^{\mathrm{2}}}}.\end{array}$$

Figure 10 Illustration of the absolute (left) and relative (right) components of the retrieval uncertainty. The top panels illustrate their dependence on thermal contrast, the bottom panels show the normalized count. Data in this plot originate from IASI-B observations on 15 January, April, July and October 2021, morning overpass, land only and between 60^∘ S and 60^∘ N. The approximations from Eqs. (72) and (76) are shown in black in the top panels.

7.4.1 Absolute uncertainty contribution

The first term is in the optically thin limit independent of the HRI and thus the column, and solely depends on the scene conditions:

$$\begin{array}{}\text{(71)}& {\mathit{\sigma }}_{\text{abs}}={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{abs},\mathrm{HRI}}}{\left|\text{SF}\right|}}={\displaystyle \frac{\sqrt{\mathrm{1}+{\mathrm{0.1}}^{\mathrm{2}}}}{\left|\text{SF}\right|}}\approx {\displaystyle \frac{\mathrm{1}}{\left|\text{SF}\right|}}.\end{array}$$

It is this term that is used as part of the post-filter to determine whether there is enough intrinsic sensitivity (thermal contrast) to make a valid measurement, i.e. one whose uncertainty is not completely overwhelmed by the instrumental noise. Currently, the post-filter threshold is set to ${\mathit{\sigma }}_{\mathrm{abs}}<\mathrm{1.5}\times {\mathrm{10}}^{\mathrm{16}}$ molec. cm⁻². Note also that a scene-dependent detection threshold of the measurements (typically taken as HRI >3) is conveniently expressed in terms of the absolute uncertainty as X_thres=3σ_abs.

The absolute uncertainty contribution is illustrated in the left panels of Fig. 10 for the IASI morning overpass (land observations between 60^∘ S and 60^∘ N) as a function of thermal contrast (TC). As before, we define TC as the brightness temperature of the surface minus the temperature at half the boundary layer height. The absolute uncertainty starts from around 1×10¹⁵ molec. cm⁻² and increases as expected with decreasing thermal contrast, with a global median of 4×10¹⁵ molec. cm⁻². Observing the inverse proportionality with thermal contrast, the following empirical formula can be used to obtain ballpark estimates of the absolute uncertainty or sensitivity of the IASI NH₃ retrieval (for positive thermal contrasts):

$$\begin{array}{}\text{(72)}& {\mathit{\sigma }}_{\mathrm{abs}}={\displaystyle \frac{\mathrm{3.6}\times {\mathrm{10}}^{\mathrm{16}}}{\text{TC}}}{\displaystyle \frac{\mathrm{molec}.\mathrm{K}}{{\mathrm{cm}}^{\mathrm{2}}}}.\end{array}$$

The constants were determined from a fit of the data shown in Fig. 10. Expressed in terms of Q20 and Q80 quantiles the estimated absolute retrieval uncertainty of IASI (mid-latitude, land, morning overpass) can also be summarized as

$$\begin{array}{}\text{(73)}& {\mathit{\sigma }}_{\mathrm{abs}}=[\mathrm{2.5}-\mathrm{6.6}]\times {\mathrm{10}}^{\mathrm{15}}\mathrm{molec}.{\mathrm{cm}}^{-\mathrm{2}}.\end{array}$$

7.4.2 Relative uncertainty contribution

The term σ_rel is proportional to the column and hence expresses a relative uncertainty for fixed atmospheric conditions,

$$\begin{array}{}\text{(74)}& {\mathit{\sigma }}_{\mathrm{rel}}=\sqrt{{\mathrm{0.1}}^{\mathrm{2}}+{\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{SF}}^{\mathrm{2}}}{{\text{SF}}^{\mathrm{2}}}}}(X-B)\approx {\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{SF}}}{\left|\text{SF}\right|}}(X-B),\end{array}$$

or

$$\begin{array}{}\text{(75)}& {\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{rel}}}{(X-B)}}={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{SF}}}{\left|\text{SF}\right|}}.\end{array}$$

This term is illustrated in the right panels of Fig. 10. Again, we observe an inverse proportionality with thermal contrast, which can be approximated as

$$\begin{array}{}\text{(76)}& {\mathit{\sigma }}_{\mathrm{rel}}=\left(\mathrm{0.07}+{\displaystyle \frac{\mathrm{1.6}\mathrm{K}}{\text{TC}}}\right){\mathrm{NH}}_{\mathrm{3}}.\end{array}$$

For typical morning land observations, the relative contribution to the uncertainty starts from around 14 % (corresponding to a TC of 20 K). Expressed in terms of Q20 and Q80 quantiles, the estimated relative retrieval uncertainty of ANNI (mid-latitude, land, morning overpass) can be summarized as

$$\begin{array}{}\text{(77)}& {\mathit{\sigma }}_{\mathrm{rel}}=[\mathrm{19}-\mathrm{36}]\mathit{\%}{\mathrm{NH}}_{\mathrm{3}}.\end{array}$$

8 Conclusions

In this paper, we presented v4 of the NH₃ ANNI retrieval. The most important change is the introduction of averaging kernels, which will greatly ease future model assimilation and comparisons with independent measurements or model output. Most other changes to ANNI v4 contribute to the overall consistency of the product. An example is the incorporation of the temporally consistent cloud flag. The improved way of calculating the HRI makes the product more robust across the different IASI instruments and more temporally harmonious. Importantly, the HRI became also less sensitive to small errors in the forward model related to the instrumental line shape function. Previous versions were biased low by some 10 %–20 % due to such errors. Theoretically we can now exclude the existence of large biases of this sort. We also demonstrate this with an optimal estimation experiment. In addition to the AVKs, we revised the uncertainty calculation and now provide better and more comprehensive information on the expected error of the measurement. We also show how the retrieval uncertainty contains a part proportional to the column and a part that is independent of the column. In the near future, the most important changes will gradually be implemented for all the other tracers retrieved with ANNI (AVKs, the use of generalized covariance matrices and the better treatment of uncertainties).