The Advanced Infra-Red WAter Vapour Estimator (AIRWAVE) version 2: algorithm evolution, dataset description and performance improvements

Total column water vapour (TCWV) is a key atmospheric variable which is generally evaluated on global scales through the use of satellite data. Recently a new algorithm, called AIRWAVE (Advanced Infra-Red WAter Vapour Estimator), has been developed for the retrieval of the TCWV from the Along-Track Scanning Radiometer (ATSR) instrument series. The AIRWAVE algorithm retrieves TCWV by exploiting the dual view of the ATSR instruments using the infrared channels at 10.8 and 12 µm and nadir and forward observation geometries. The algorithm was used to produce a TCWV database over sea from the whole ATSR mission. When compared to independent TCWV products, the AIRWAVE version 1 (AIRWAVEv1) database shows very good agreement with an overall bias of 3 % all over the ATSR missions. A large contribution to this bias comes from the polar and the coastal regions, where AIRWAVE underestimates the TCWV amount. In this paper we describe an updated version of the algorithm, specifically developed to reduce the bias in these regions. The AIRWAVE version 2 (AIRWAVEv2) accounts for the atmospheric variability at different latitudes and the associated seasonality. In addition, the dependency of the retrieval parameters on satellite across-track viewing angles is now explicitly handled. With the new algorithm we produced a second version of the AIRWAVE dataset. As for AIRWAVEv1, the quality of the AIRWAVEv2 dataset is assessed through the comparison with the Special Sensor Microwave/Imager (SSM/I) and with the Analyzed RadioSounding Archive (ARSA) TCWV data. Results show significant improvements in both biases (from 0.72 to 0.02 kg m⁻²) and standard deviations (from 5.75 to 4.69 kg m⁻²), especially in polar and coastal regions. A qualitative and quantitative estimate of the main error sources affecting the AIRWAVEv2 TCWV dataset is also given. The new dataset has also been used to estimate the water vapour climatology from the 1991–2012 time series.

1 Introduction

A key issue in assessing the climate change is the precise knowledge of the distribution and variability of the total column of water vapour (TCWV), i.e. the vertically integrated atmospheric water vapour content. Actually, TCWV is closely linked to clouds, precipitation and thus to the hydrological cycle (Allan et al., 2014). For this reason it is one of the GCOS (Global Climate Observing System) essential climate variables (ECVs). Since water vapour plays such a crucial role in meteorological as well as in climatological aspects, it is important to gather spatial and temporal thorough information about its distribution. On a global scale, this can be achieved through the use of satellite missions. In the last decades measurements from several sensors were used for this purpose. Among them, sensors operating in the microwave region, such as the Special Sensor Microwave Imager (SSM/I on board Defense Meteorological Satellite Program, DMSP, satellites), are used to infer the accurate amount of TCWV over ocean surfaces (Wentz, 1997), while sensors operating in the visible and near-infrared spectral range provide precise TCWV retrieval on land surfaces (e.g. the ENVISAT Medium Resolution Imaging Spectrometer (MERIS), Lindstrot et al., 2012, or Moderate Resolution Imaging Spectroradiometer, MODIS, on the Terra and Aqua satellites, Diedrich et al., 2015).

TCWV retrievals from infrared spectral regions were performed from Advanced Very High Resolution Radiometer (AVHRR; Emery, 1992) measurements, using the split window technique, (Li et al., 2003; Sobrino et al., 1991) and from MODIS (Seemann et al., 2003). TCWV retrievals from infrared channels over land suffer from the limited knowledge of the temperature and the emissivity of land surfaces (Lindstrot et al., 2014). The Along-Track Scanning Radiometer (ATSR, Delderfield et al., 1986) instrument series had as its main objective the accurate retrieval of sea surface temperature for climate studies. However, Casadio et al. (2016) demonstrated that it is possible to retrieve accurate and precise TCWV from its daytime and night-time measurements, using the ATSR brightness temperature (BT) collected from nadir and forward views in the channels at 10.8 and 12 µm in clear-sky daytime and night-time sea scenes. The algorithm (named AIRWAVE, Advanced Infra-Red WAter Vapour Estimator) exploits a sea emissivity dataset and calculations made with a dedicated radiative transfer model (RTM). A detailed description of the AIRWAVE algorithm is given in Casadio et al. (2016). The first version of the AIRWAVE TCWV dataset (hereafter AIRWAVEv1), spanning from 1991 to 2012, is freely available from the GEWEX G-VAP website (G-VAP, 2018) in the form of monthly fields at $\mathrm{2}{}^{\circ }\times \mathrm{2}{}^{\circ }$ regular grid resolution from 2003 to 2008 (Schroder et al., 2018). Due to the legacy of the ATSR series, and the fact that the radiances are a fundamental climate dataset record, the AIRWAVE dataset is an important resource for water vapour studies. It is worth underlying here that AIRWAVEv1 was developed to demonstrate the possibility of retrieving TCWV values from the ATSR measurements. The main goal pursued in its development was to have simple software that could produce good results when compared to independent datasets. For this reason, in the AIRWAVEv1 algorithm several approximations were made. AIRWAVEv1 use fixed retrieval parameters along the globe and TCWV are corrected for viewing angle variability in nadir and slant by using an empirical correction factor.

Papandrea et al. (2018), aiming to validate the AIRWAVEv1 dataset, compared the data with the TCWV from SSM/I and Analyzed RadioSounding Archive, ARSA (2018) for the whole mission. This exercise demonstrated that the AIRWAVEv1 dataset was in general of high quality (average correlative bias of 0.72 kg m⁻² vs. SSM/I and 0.80 kg m⁻² vs. ARSA, below the 1 kg m⁻² indicated in the GlobVapour project, Lindstrot et al., 2010) apart for the polar regions and some coastal regions, where an underestimation of the TCWV was found. In this paper we describe the new version of the AIRWAVE algorithm (hereafter AIRWAVEv2) developed to overcome these weaknesses by accounting for latitudinal and angular variations of the retrieval parameters. The new algorithm has been applied to all the available ATSR level 1B top-of-atmosphere radiance products acquired over water surfaces in clear-sky and in day and night conditions (same as for AIRWAVEv1) to produce the AIRWAVEv2 dataset. Here, we show the new TCWV climatologies derived from 20 years of ATSR data together with the results of an extensive validation exercise repeating the same comparisons reported in Papandrea et al. (2018). The new dataset shows improvements in terms of both bias and spread of the differences with respect to another dataset and what was achieved with AIRWAVEv1.

This article is structured as follows: in Sect. 2 we describe the new algorithm developed to produce AIRWAVEv2, the improvements in retrieval scenarios and the strategy used for the selection of latitude- and seasonal-dependent retrieval parameters. In Sect. 3 we describe the AIRWAVEv2 dataset, the TCWV climatology and its validation against SSM/I and ARSA data and compare the performances of AIRWAVEv2 and AIRWAVEv1. Finally, conclusions are given in Sect. 4.

2 AIRWAVE version 2

Papandrea et al. (2018) demonstrated the high quality of AIRWAVEv1 by comparing the retrieved TCWV with corresponding SSM/I and ARSA TCWV. However, in the same paper, the authors highlighted that, at latitudes higher than 50^∘, the agreement was not as good as for the rest of the globe. They speculated that this was due to the fact that AIRWAVEv1 makes use of retrieval parameters calculated though RTM simulations of tropical and midlatitude atmospheric scenarios, which are then averaged and used for the whole globe. This choice was driven by the consideration that, AIRWAVEv1 being applicable to water and cloud-free scenes only, the number of cloud-free measurements over the sea at high latitudes is significantly smaller than at midlatitudes and tropical regions. Thus, a trade-off between generality (i.e. good precision at all latitudes), actual latitudinal coverage of cloud-free measurements and software complexity was the main driver for this choice. Moreover, AIRWAVEv1 makes use of the retrieval parameters computed for the along-track viewing geometries only and uses an a posteriori correction for the scenes pointing outside the orbit track.

The need to have a TCWV dataset of homogeneous quality at all latitudes and viewing geometries has driven the development of an improved version of the AIRWAVE algorithm, AIRWAVEv2. The improvements were achieved through three main steps. Firstly, we modified the way in which some of the approximations of the solving equations were handled, leading to an improved retrieval precision. Secondly, we compute the retrieval parameters for different latitude bands and for 4 months that, in the retrieval, are used as look-up-tables. Finally, we calculated the retrieval parameters for different viewing angles to directly account for across-track variations. We discuss these modifications in the following subsections. All the computations described in the paper have been made using the HITRAN2008 database (Rothman et al., 2009) for the spectroscopic data and the IG2 database (Remedios et al., 2007) version 4.1 for the atmospheric scenarios. The IG2 database was developed to be used as a model atmosphere in the analysis of the ENVISAT Michelson Interferometer for Passive Atmospheric Sounding (MIPAS) measurements (that cover the same spectral region of the IR channels of ATSR). The IG2 database contains atmospheric vertical profiles of pressure, temperature and abundances of the molecules active in the MIPAS spectral region, which are different for each year of the mission and divided into six latitudinal bands (polar, midlatitude and equatorial for the northern and southern hemispheres) and four seasons.

2.1 Improvements in the solving equations

The starting point for the calculations of AIRWAVEv2 retrieval parameters is the master equation of AIRWAVE version 1. Since the expressions are the same for both nadir and forward geometry, here we report the equations for the general case, omitting geometries denoted by subscript NAD (for NADIR) and FWD (for FORWARD). Here we reproduce some of the equations reported in Casadio et al. (2016) to help the reader in the comprehension of this article. The master equation of the AIRWAVE algorithm is Eq. (12) of the above-mentioned work (now Eq. 1):

$$\begin{array}{}\text{(1)}& {\displaystyle }\ln {\displaystyle \frac{J_{\mathrm{1}}^{{\mathit{\lambda }}_{\mathrm{1}}}}{J_{\mathrm{2}}^{{\mathit{\lambda }}_{\mathrm{2}}}}}=\ln {\displaystyle \frac{F_{\mathrm{1}}^{{\mathit{\lambda }}_{\mathrm{1}}}}{F_{\mathrm{2}}^{{\mathit{\lambda }}_{\mathrm{2}}}}}+\ln {\displaystyle \frac{{\mathit{ϵ}}_{\mathrm{1}}^{{\mathit{\lambda }}_{\mathrm{1}}}}{{\mathit{ϵ}}_{\mathrm{2}}^{{\mathit{\lambda }}_{\mathrm{2}}}}}+\ln {\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{1}}^{{\mathit{\lambda }}_{\mathrm{1}}}}{{\mathit{\gamma }}_{\mathrm{2}}^{{\mathit{\lambda }}_{\mathrm{2}}}}}+{\mathit{\lambda }}_{\mathrm{2}}{\mathit{\tau }}_{\mathrm{2}}-{\mathit{\lambda }}_{\mathrm{1}}{\mathit{\tau }}_{\mathrm{1}}.\end{array}$$

The subscript numbers 1 and 2 represent the terms calculated in the 10.8 and 12 µm channels respectively. λ₁ is the value of the frequency in the 10.8 µm channel and λ₂ is that in the 12 µm channel. J₁ is the radiance that reaches the TOA for the 10.8 µm channel, J₂ is the radiance that reaches the TOA for the 12 µm channel, F includes the atmospheric (J_a) and surface radiance (J_s) contribution and is $F=\mathrm{1}+\frac{J_{\mathrm{a}}}{e^{-\mathit{\tau }}{J}_{\mathrm{s}}}$, ϵ is the sea emissivity, γ is a constant arising from the Planck law, τ are the optical depths at the two wavelengths. Since only H₂O and CO₂ significantly affect the optical depth into ATSR thermal infrared (TIR) channels we can write the following:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\ln {\displaystyle \frac{J_{\mathrm{1}}^{{\mathit{\lambda }}_{\mathrm{1}}}}{J_{\mathrm{2}}^{{\mathit{\lambda }}_{\mathrm{2}}}}}=\ln {\displaystyle \frac{F_{\mathrm{1}}^{{\mathit{\lambda }}_{\mathrm{1}}}}{F_{\mathrm{2}}^{{\mathit{\lambda }}_{\mathrm{2}}}}}+\ln {\displaystyle \frac{{\mathit{ϵ}}_{\mathrm{1}}^{{\mathit{\lambda }}_{\mathrm{1}}}}{{\mathit{ϵ}}_{\mathrm{2}}^{{\mathit{\lambda }}_{\mathrm{2}}}}}+\ln {\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{1}}^{{\mathit{\lambda }}_{\mathrm{1}}}}{{\mathit{\gamma }}_{\mathrm{2}}^{{\mathit{\lambda }}_{\mathrm{2}}}}}+{\mathit{\lambda }}_{\mathrm{2}}{\mathit{\tau }}_{\mathrm{2}}^{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}\\ \text{(2)}& {\displaystyle }& {\displaystyle }-{\mathit{\lambda }}_{\mathrm{1}}{\mathit{\tau }}_{\mathrm{1}}^{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}+{\mathit{\lambda }}_{\mathrm{2}}{\mathit{\tau }}_{\mathrm{2}}^{{\mathrm{CO}}_{\mathrm{2}}}-{\mathit{\lambda }}_{\mathrm{1}}{\mathit{\tau }}_{\mathrm{1}}^{{\mathrm{CO}}_{\mathrm{2}}}.\end{array}$$

By renaming,

$$\begin{array}{}\text{(3)}& {\displaystyle }G=\ln {\displaystyle \frac{F_{\mathrm{1}}^{{\mathit{\lambda }}_{\mathrm{1}}}}{F_{\mathrm{2}}^{{\mathit{\lambda }}_{\mathrm{2}}}}},E=\ln {\displaystyle \frac{{\mathit{ϵ}}_{\mathrm{1}}^{{\mathit{\lambda }}_{\mathrm{1}}}}{{\mathit{ϵ}}_{\mathrm{2}}^{{\mathit{\lambda }}_{\mathrm{2}}}}},\mathit{\chi }=\ln {\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{1}}^{{\mathit{\lambda }}_{\mathrm{1}}}}{{\mathit{\gamma }}_{\mathrm{2}}^{{\mathit{\lambda }}_{\mathrm{2}}}}},\end{array}$$

we get

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\ln {\displaystyle \frac{J_{\mathrm{1}}^{{\mathit{\lambda }}_{\mathrm{1}}}}{J_{\mathrm{2}}^{{\mathit{\lambda }}_{\mathrm{2}}}}}=G+E+\mathit{\chi }+{\mathit{\lambda }}_{\mathrm{2}}{\mathit{\tau }}_{\mathrm{2}}^{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}-{\mathit{\lambda }}_{\mathrm{1}}{\mathit{\tau }}_{\mathrm{1}}^{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}\\ \text{(4)}& {\displaystyle }& {\displaystyle }+{\mathit{\lambda }}_{\mathrm{2}}{\mathit{\tau }}_{\mathrm{2}}^{{\mathrm{CO}}_{\mathrm{2}}}-{\mathit{\lambda }}_{\mathrm{1}}{\mathit{\tau }}_{\mathrm{1}}^{{\mathrm{CO}}_{\mathrm{2}}}.\end{array}$$

In AIRWAVEv1 we assume that the optical depth (τ) is the product of the vertical column of H₂O using the relative effective absorption cross section (λσ), normalised to the air mass factor (AMF) for the given line-of-sight angle:

$$\begin{array}{}\text{(5)}& {\displaystyle }{\mathit{\lambda }}_{\mathrm{2}}{\mathit{\tau }}_{\mathrm{2}}^{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}-{\mathit{\lambda }}_{\mathrm{1}}{\mathit{\tau }}_{\mathrm{1}}^{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}={\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{2}}-{\mathit{\lambda }}_{\mathrm{1}}{\mathit{\sigma }}_{\mathrm{1}}}{\mathrm{AMF}}}\mathrm{TCWV}.\end{array}$$

This equation shows that linear behaviour exists between the water vapour optical depth and the TCWV. The linear dependence is exploited to solve the AIRWAVE equation and to retrieve TCWV.

In the development of AIRWAVEv2 we investigated the possibility of finding a more accurate solution to the AIRWAVE equation while still preserving the linear dependence between water vapour optical depth and TCWV. We recall here that in AIRWAVEv1 the water absorption cross sections were obtained using MODTRAN (Berk et al., 2008), while all the other values were obtained with the dedicated RTM, developed for ATSR measurements simulations and described in Casadio et al. (2016). For AIRWAVEv2, a different approach was adopted for the calculation of effective absorption cross section. We have simulated ATSR synthetic radiances for the different atmospheric scenarios of the IG2 database and thus with different water vapour contents. A detailed description of these simulations is given in Sect. 2.2.

Using these simulations and ATSR-SSM/I-collocated TCWV, we verified that $\ln \frac{J_{\mathrm{1}}^{{\mathit{\lambda }}_{\mathrm{1}}}}{J_{\mathrm{2}}^{{\mathit{\lambda }}_{\mathrm{2}}}}$ correctly reproduces the real measurement behaviour as a function of TCWV and that this relation is a first linear approximation.

In Fig. 2 the coloured dots represent the values of the logarithm of the radiance ratio in Eq. (1) as a function of the TCWV for the different atmospheric scenarios. We report only the values obtained for the subsatellite scans using the IG2 water profiles for the summer season multiplied for 0.5 and 1.5. The different colours represent different latitude bands (going from red for tropical to blue for polar). The grey dots represent the radiance ratio calculated from along-track AATSR measurements on the 5 and 6 August 2008 aggregated at SSM/I resolution (0.25${}^{\circ }\times \mathrm{0.25}{}^{\circ }$). The value of TCWV associated with each AATSR subsatellite measurement was obtained from coincident SSM/I measurements. In order to minimise the impact of random error, only measurements with SSM/I pixel coverage (calculated as the ratio between the actual and the maximum number of ATSR measurements that can be present into a SSM/I pixel) greater than 10 % were used for this exercise. Figure 2 shows that (a) the simulated radiances correctly reproduce the real measurement behaviour, and (b) the relation between the radiances and the TCWV can be considered linear. Actually, in this case we find a correlation of 0.904 for real data and 0.92 for the simulated ones (p value of $\mathrm{7.3}\times {\mathrm{10}}^{-\mathrm{05}}$).

Therefore we can now rewrite Eq. (4), isolating the terms that account for the water content (${\mathit{\tau }}_{\mathrm{1}}^{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ and ${\mathit{\tau }}_{\mathrm{2}}^{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$):

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\ln {\displaystyle \frac{J_{\mathrm{1}}^{{\mathit{\lambda }}_{\mathrm{1}}}}{J_{\mathrm{2}}^{{\mathit{\lambda }}_{\mathrm{2}}}}}-G-\mathit{\chi }-E-{\mathit{\lambda }}_{\mathrm{2}}{\mathit{\tau }}_{\mathrm{2}}^{{\mathrm{CO}}_{\mathrm{2}}}+{\mathit{\lambda }}_{\mathrm{1}}{\mathit{\tau }}_{\mathrm{1}}^{{\mathrm{CO}}_{\mathrm{2}}}\\ \text{(6)}& {\displaystyle }& {\displaystyle }={\mathit{\lambda }}_{\mathrm{2}}{\mathit{\tau }}_{\mathrm{2}}^{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}-{\mathit{\lambda }}_{\mathrm{1}}{\mathit{\tau }}_{\mathrm{1}}^{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}=\mathrm{\Delta }\mathit{\tau }.\end{array}$$

The G term of Eqs. (3) and (4) is not as constant as supposed, and as is partially verified in Casadio et al. (2016) and may depend on the differing water vapour content. For this reason, for each atmospheric scenario the average of all the G values obtained with different water vapour contents is used (to vary the water vapour content we multiplied the water vapour profile for 0.5, 0.75, 1., 1.25 and 1.5). Equation (6) can thus be written as follows:

$$\begin{array}{}\text{(7)}& {\displaystyle }\ln {\displaystyle \frac{J_{\mathrm{1}}^{{\mathit{\lambda }}_{\mathrm{1}}}}{J_{\mathrm{2}}^{{\mathit{\lambda }}_{\mathrm{2}}}}}-G_{\mathrm{AVG}}-\mathit{\chi }-E-{\mathit{\lambda }}_{\mathrm{2}}{\mathit{\tau }}_{\mathrm{2}}^{{\mathrm{CO}}_{\mathrm{2}}}+{\mathit{\lambda }}_{\mathrm{1}}{\mathit{\tau }}_{\mathrm{1}}^{{\mathrm{CO}}_{\mathrm{2}}}=\mathrm{\Delta }\mathit{\tau }.\end{array}$$

Therefore for each scenario and geometry we can write

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\ln {\displaystyle \frac{J_{\mathrm{1}}^{{\mathit{\lambda }}_{\mathrm{1}}}}{J_{\mathrm{2}}^{{\mathit{\lambda }}_{\mathrm{2}}}}}-G_{\mathrm{AVG}}-\mathit{\chi }-E-{\mathit{\lambda }}_{\mathrm{2}}{\mathit{\tau }}_{\mathrm{2}}^{{\mathrm{CO}}_{\mathrm{2}}}+{\mathit{\lambda }}_{\mathrm{1}}{\mathit{\tau }}_{\mathrm{1}}^{{\mathrm{CO}}_{\mathrm{2}}}\\ \text{(8)}& {\displaystyle }& {\displaystyle }=\mathrm{\Delta }\mathit{\tau }=\mathrm{\Delta }\mathit{\sigma }\cdot \mathrm{TCWV}+\mathrm{\Delta }\mathit{\rho }.\end{array}$$

In this equation, Δσ and Δρ represent the slope and intercept of the straight line representing the behaviour of the term containing the radiances as a function of the TCWV. In the testing version of the AIRWAVEv2 code, we estimated these parameters using the values of the radiances and the TCWV obtained perturbing the IG2 water vapour amount by factors of 0.5 and 1.5. By grouping the terms in Eq. (8) as in Casadio et al. (2016) we get

$$\begin{array}{}\text{(9)}& {\displaystyle }\mathrm{TCWV}=\mathrm{\Phi }-{\displaystyle \frac{G}{\mathrm{\Delta }\mathit{\sigma }}},\end{array}$$

where Φ is the “water vapour pseudo-column” that in AIRWAVEv2 is defined as follows:

$$\begin{array}{}\text{(10)}& {\displaystyle }\mathrm{\Phi }={\displaystyle \frac{\ln \frac{J_{\mathrm{1}}^{{\mathit{\lambda }}_{\mathrm{1}}}}{J_{\mathrm{2}}^{{\mathit{\lambda }}_{\mathrm{2}}}}-\mathit{\chi }-E-{\mathit{\lambda }}_{\mathrm{2}}{\mathit{\tau }}_{\mathrm{2}}^{{\mathrm{CO}}_{\mathrm{2}}}+{\mathit{\lambda }}_{\mathrm{1}}{\mathit{\tau }}_{\mathrm{1}}^{{\mathrm{CO}}_{\mathrm{2}}}-\mathrm{\Delta }\mathit{\rho }}{\mathrm{\Delta }\mathit{\sigma }}}.\end{array}$$

This formula is slightly different from the one used in Casadio et al. (2016) due to the presence of the Δρ term:

$$\begin{array}{}\text{(11)}& {\displaystyle }\mathrm{\Phi }={\displaystyle \frac{\ln \frac{J_{\mathrm{1}}^{{\mathit{\lambda }}_{\mathrm{1}}}}{J_{\mathrm{2}}^{{\mathit{\lambda }}_{\mathrm{2}}}}-\mathit{\chi }-E-{\mathit{\lambda }}_{\mathrm{2}}{\mathit{\tau }}_{\mathrm{2}}^{{\mathrm{CO}}_{\mathrm{2}}}+{\mathit{\lambda }}_{\mathrm{1}}{\mathit{\tau }}_{\mathrm{1}}^{{\mathrm{CO}}_{\mathrm{2}}}}{\mathrm{\Delta }\mathit{\sigma }}}.\end{array}$$

If in Eq. (9) now we make explicit the dependence on the viewing angles, we get the following:

$$\begin{array}{}\text{(12)}& {\displaystyle }& {\displaystyle }{\mathrm{TCWV}}_{\mathrm{NAD}}={\mathrm{\Phi }}_{\mathrm{NAD}}-{\displaystyle \frac{G_{\mathrm{NAD}}}{\mathrm{\Delta }{\mathit{\sigma }}_{\mathrm{NAD}}}}\text{(13)}& {\displaystyle }& {\displaystyle }\mathrm{and}{\mathrm{TCWV}}_{\mathrm{FWD}}={\mathrm{\Phi }}_{\mathrm{FWD}}-{\displaystyle \frac{G_{\mathrm{FWD}}}{\mathrm{\Delta }{\mathit{\sigma }}_{\mathrm{FWD}}}},\end{array}$$

with

$$\begin{array}{}\text{(14)}& {\displaystyle }& {\displaystyle }{\mathrm{\Phi }}_{\mathrm{NAD}}={\displaystyle \frac{\left\{\begin{array}{c}\ln \frac{J_{{\mathrm{1}}_{\mathrm{NAD}}}^{{\mathit{\lambda }}_{\mathrm{1}}}}{J_{{\mathrm{2}}_{\mathrm{NAD}}}^{{\mathit{\lambda }}_{\mathrm{2}}}}-\mathit{\chi }-E_{\mathrm{NAD}}-{\mathit{\lambda }}_{\mathrm{2}}{\mathit{\tau }}_{{\mathrm{2}}_{\mathrm{NAD}}}^{{\mathrm{CO}}_{\mathrm{2}}}\\ +{\mathit{\lambda }}_{\mathrm{1}}{\mathit{\tau }}_{{\mathrm{1}}_{\mathrm{NAD}}}^{{\mathrm{CO}}_{\mathrm{2}}}-\mathrm{\Delta }{\mathit{\rho }}_{\mathrm{NAD}}\end{array}\right\}}{\mathrm{\Delta }{\mathit{\sigma }}_{\mathrm{NAD}}}}\text{(15)}& {\displaystyle }& {\displaystyle }{\mathrm{\Phi }}_{\mathrm{FWD}}={\displaystyle \frac{\left\{\begin{array}{c}\ln \frac{J_{{\mathrm{1}}_{\mathrm{FWD}}}^{{\mathit{\lambda }}_{\mathrm{1}}}}{J_{{\mathrm{2}}_{\mathrm{FWD}}}^{{\mathit{\lambda }}_{\mathrm{2}}}}-\mathit{\chi }-E_{\mathrm{FWD}}-{\mathit{\lambda }}_{\mathrm{2}}{\mathit{\tau }}_{{\mathrm{2}}_{\mathrm{FWD}}}^{{\mathrm{CO}}_{\mathrm{2}}}\\ +{\mathit{\lambda }}_{\mathrm{1}}{\mathit{\tau }}_{{\mathrm{1}}_{\mathrm{FWD}}}^{{\mathrm{CO}}_{\mathrm{2}}}-\mathrm{\Delta }{\mathit{\rho }}_{\mathrm{FWD}}\end{array}\right\}}{\mathrm{\Delta }{\mathit{\sigma }}_{\mathrm{FWD}}}}.\end{array}$$

The TCWV can be estimated through knowledge of G. Actually, for single-view geometry, the presence of the G term can affect the accurate determination of TCWV. Thus, the G variability suggests that it will be desirable to avoid this term in TCWV derivation. In AIRWAVE this is done by exploiting the dual-view capability of the ATSR instruments and by assuming a perfect collocation between NAD and FWD measurements:

$$\begin{array}{}\text{(16)}& {\displaystyle }\mathrm{TCWV}=\mathit{\alpha }\cdot {\mathrm{\Phi }}_{\mathrm{NAD}}+\mathit{\beta }\cdot {\mathrm{\Phi }}_{\mathrm{FWD}},\end{array}$$

where

$$\begin{array}{}\text{(17)}& {\displaystyle }& {\displaystyle }\mathit{\alpha }={\displaystyle \frac{\mathrm{1}}{\mathrm{1}-\frac{\mathrm{\Delta }{\mathit{\sigma }}_{\mathrm{FWD}}}{\mathit{\delta }\cdot \mathrm{\Delta }{\mathit{\sigma }}_{\mathrm{NAD}}}}}\mathrm{and}\mathit{\beta }={\displaystyle \frac{\mathrm{1}}{\mathrm{1}-\frac{\mathit{\delta }\cdot \mathrm{\Delta }{\mathit{\sigma }}_{\mathrm{NAD}}}{\mathrm{\Delta }{\mathit{\sigma }}_{\mathrm{FWD}}}}}\text{(18)}& {\displaystyle }& {\displaystyle }\mathrm{with}\mathit{\delta }\approx {\displaystyle \frac{G_{\mathrm{FWD}}}{G_{\mathrm{NAD}}}}.\end{array}$$

Equations (10), (12) and (16) are the solving equations used in AIRWAVEv2, while AIRWAVEv1 makes use of Eqs. (11), (12) and (16). Equations  (10), (12) and (16) were solved for the 11 couples of viewing angles corresponding to the tie points. The angles cover a range from 0 to 21^∘ in the NAD case and from 53 to 55^∘ in the FWD case. The new equations were used to compute a new set of retrieval parameters. For consistency purposes, in AIRWAVEv2 they were computed with the dedicated RTM as they were for AIRWAVEv1. An example of the difference between the parameters used for AIRWAVEv1 and AIRWAVEv2 is given in Table 1 for the tropical scenario and subsatellite view configuration. As can be noticed, the larger differences are for α, β and Δσ parameters: there is a reduction in the α and β parameters by a factor of about 50 from AIRWAVEv1 to AIRWAVEv2, while Δσ is reduced by a factor of 30.

Table 1 AIRWAVE v1 and AIRWAVEv2 retrieval parameters for tropical summer atmosphere and along-track configuration (subsatellite view at 0 and 55^∘).

These changes have a direct effect on the retrieval precision. We can have an estimate of the improvements of the precision of AIRWAVEv2 retrievals against AIRWAVEv1 using the following consideration: in AIRWAVE we can estimate the expected precision by multiplying the measurement random error by a factor of α∕Δσ. Therefore, in AIRWAVEv1, we multiply the random error by a factor of $\mathrm{53}/\mathrm{1.9}(=\mathrm{28})$, while in AIRWAVEv2 the multiplicative factor is $\mathrm{1.7}/\mathrm{0.08}(=\mathrm{21})$. In the case of polar atmosphere this ratio is further reduced, reaching a value of 13. Since the overall random error is about 0.25 % for AATSR and ATSR-2 and 0.6 % for ATSR-1, the precision for AATSR and ATSR-2 improves from 7 % to 5 % and for ATRS-1 from 17 % to 12 % for tropical atmosphere. In AIRWAVEv2 it reaches 3 % for AATSR and ATSR-2 and 8 % for ATSR-1 in polar atmosphere. In the worst cases, the precision has at maximum a 1 % higher value in the extreme across-track of the swath with respect to the subsatellite points. From these considerations it follows that AIRWAVEv2 parameters should improve the retrieval.

2.2 Improvements in the retrieval scenario

In AIRWAVEv1 we make use of the same set of retrieval parameters for all measurements of a single ATSR instrument. Independent sets of parameters are calculated for the three missions, while within each mission no dependencies on different atmospheric–surface conditions or seasons was considered. In AIRWAVEv2 the retrieval parameters are estimated, not only according to the instrument type but also accounting for possible latitudinal and seasonal variations. To this aim, we have used the aforementioned RTM to compute all the required quantities by exploiting the model atmospheres of the IG2 database. Since the IG2 database was specifically developed for the ENVISAT MIPAS mission, it covers only the time range from 2002 to 2012, while the ATSR series operated from 1991 to 2012. However, the interannual variations of most of the species active in the ATSR thermal infrared spectral range are generally much smaller than the corresponding seasonal ones. Therefore we used the data for 1 year only (2010), and we considered the interannual variations to be systematic error sources (see Appendix A for an estimate of these errors).

To better reproduce the variability of the atmospheric scenarios that were observed by the ATSR instruments, we therefore calculated the retrieval parameters exploiting the profiles for all six latitude bands and the four seasons included in the IG2 datasets for 1 year only.

The AIRWAVE solving equations do not explicitly make use of the sea surface temperature (SST). However, SST is used in the radiative transfer (RT) computations to estimate the retrieval parameters. The sea surface emissivity values are instead used both in the parameter estimation and in the AIRWAVE retrieval. While over land the emissivity is characterised by a large spatial difference (it indeed varies as a function of soil type, vegetation cover, etc.), over sea its variation is in general relatively small. For this reason, in AIRWAVEv1, we used constant emissivity values calculated for the nadir and forward viewing angles with fixed SST (285 K) and wind speed (3 m s⁻¹). In the new version of the algorithm, coherently with the approach used for the atmospheric scenarios, the retrieval parameters have been computed using dedicated SST values for each season and latitude band. The SST monthly means were produced for the corresponding six latitude bands and for the four seasons using ECMWF ERA-Interim daily field data with a regular latitude–longitude grid of $\mathrm{0.75}{}^{\circ }\times \mathrm{0.75}{}^{\circ }$ (241×480 grid points).

The emissivity of each scenario has been computed using data extracted from the University of Edinburgh database (Embury et al., 2008). This dataset contains emissivities tabulated as a function of wave number (600–3350 cm⁻¹ or 3–16.7 µm), viewing angle (0–85^∘), temperature (270–310 K), and wind speed (0–25 m s⁻¹ at 12.5 m). For the RTM computations we used the full spectral dependency of the emissivity. Since in Eq. (3) of AIRWAVE we use a single emissivity value for each channel, we estimated it by convolving the spectral emissivity with the ATSR filter functions. The nadir and forward viewing angles of the instruments have been defined at 11 tie points of the ATSR swath (pixels associated with specific points equally spaced across a single image or instrument scan). For each tie point we then used the corresponding viewing angles to extract the correct emissivity values, with a fixed wind speed (3 m s⁻¹), as for AIRWAVEv1. Due to wind speed variability, in both algorithm versions, we prefer to fix the value of the wind speed and then treat wind variations as an error source (see Sect. A3 in Appendix).

2.3 Across-track variations of the retrieval parameters

A simplification present in AIRWAVEv1 is that the retrieval parameters are calculated only for the subsatellite viewing angles (55^∘ for forward viewing and 0^∘ as nadir viewing angle). However, due to the ATSR configuration, the nadir viewing angles vary from 0^∘ (subsatellite) to approximately 21^∘ (across-track edge of the ATSR swath, ±250 km from nadir), while the forward viewing angles range from 53 to 55^∘. Significant TCWV differences between the centre and edge swaths are thus expected.

In AIRWAVEv1 the across-track dependence of TCWV is corrected a posteriori. The correction was calculated on the basis of TCWV retrievals performed over simulated brightness temperatures (BTs).

Figure 1 shows the absolute difference between TCWV from across-track pixels and subsatellite pixels calculated for AATSR in different atmospheric scenarios (tropical and midlatitude) from synthetic measurements. To produce these retrievals, we simulated ATSR radiances at the 11 tie points for the nadir and forward views while also considering the surface emissivity variations with viewing angles. The atmospheric scenario, the TCWV and the SST were kept constant and they were exactly the same as that used for the subsatellite track case. Then we retrieved the TCWV for each of the 11 couples of BTs and computed the difference with respect to the value obtained at the subsatellite track position (which coincides with the TCWV reference value). The dotted line in Fig. 1 mimics the correction term adopted within the AIRWAVEv1 algorithm. As can be seen, the dotted line reproduces the general behaviour of the across-track TCWV dependence well. This is confirmed by the comparison made between AIRWAVEv1 TCWV and SSM/I or ECMWF TCWV (Casadio et al., 2016): on average we did not find any across-track bias, thus the general validity of this correction is confirmed.

Figure 1 Differences between TCWV calculated from across pixels and TCWV calculated from the subsatellite track as a function of across-track position calculated for different atmospheric scenarios and for AATSR (symbols). Grey dashed line corresponds to AIRWAVEv1 parameterisation.

However, as can be seen in Fig. 1, this approximation might not be sufficiently adequate depending on the used atmospheric scenario. Furthermore, a slight asymmetry with respect to the subsatellite track position is expected as the ATSR instruments are tilted about 4^∘ respect to the flying direction of the satellite. Therefore the a posteriori correction of AIRWAVEv1 cannot fully reproduce all these features. In AIRWAVEv2 we replaced the a posteriori correction: we calculated the retrieval parameters for each of the above-described tie points of the nadir and forward swaths, then we obtained the parameters at the exact ground pixel position interpolating these values and using the ground pixel across-track position.

Figure 2 Logarithm of radiance ratio at nadir as a function of TCWV in simulated atmospheric scenarios (purple and blue for polar, light and dark green for midlatitude and light and dark red for equatorial). Grey dots represent the same quantity using real AATSR radiances and coincident SSM/I TCWV for the along-track measurements on 5 and 6 August 2008.

2.4 Selection of the retrieval parameters

The computation of the retrieval parameters, described in the previous sections, produced a set of 1584 retrieval parameters for each ATSR mission (6 coefficients ×6 latitude bands ×4 seasons ×11 tie points) stored in dedicated look-up tables. In order to select the most suitable set of parameters for each ATSR measurement, a multivariate interpolation on a 3-dimensional grid (trilinear) has been applied to the six retrieval coefficients (Δσ, G, Δρ) and the emissivity for both the FWD and the NAD geometries.

3 AIRWAVEv2 dataset: description, climatology and validation

The AIRWAVEv2 TCWV data are produced in the same way as the AIRWAVEv1 dataset, by processing level 1B measurements acquired over water surfaces (sea and lakes) and in clear-sky conditions in both nadir and forward views (according to the level 1B cloud mask). The output files are saved in Interactive Data Language (IDL) binary files (.sav extension); however they can be easily converted to other formats (e.g. netcdf) upon request. The parameters contained into the files are structured in two groups: in the first one (named HIRES) the parameters are given at native resolution (1×1 km²), while in the second one (named SSM/I) the parameters are aggregated to SSM/I resolution ($\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ grid). Both groups contain the TCWV, the latitude, the longitude, the across-track index value (0–512) and a day–night flag. The SSM/I group, in addition, contains the value of the number of elements aggregated within the SSM/I grid cell and the standard deviation of the TCWV value associated with each cell.

The climatologies have been derived using all the available years and sensors of the ATSR family. Using the AIRWAVEv2 products aggregated on the $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ grid, we obtained for each month day and night TCWV averages and standard deviations. The final monthly files are also saved as IDL binary files but can be converted to other formats and are available upon request.

Figure 3 Climatology of daytime TCWV from the AIRWAVEv2 dataset for January (a), April (c), July (e) and October (g) from 1991 to 2012, and standard deviations for the same months (b, d, f, h).

Figures 3 and 4 show sketches of the climatology for January, April, July and October obtained from 20 years of ATSR daytime measurements. Similar results are obtained for night-time retrievals (not shown here). In Fig. 3 we report in panels (a), (c), (e) and (g) the TCWV global distribution at 0.25^∘ grid resolution and its standard deviation in panels (b), (d), (f) and (h). In Fig. 4 panels (a), (c), (e) and (g) we show the TCWV meridional mean together with its standard deviation, while on panels (b), (d), (f) and (h) we report the TCWV zonal mean.

Figure 4 TCWV distribution as function of latitude (median and standard deviation) (a, c, e, g) and as function of longitude (b, d, f, h).

The geographical distribution of the median values of the TCWV reflects the behaviour of general atmospheric circulation. Higher TCWV values, associated with strong convection in the Intertropical Convergence Zone (ITCZ), are located around the equator, while lower values are found in the polar regions. Also, the zonal means as a function of latitude reflects this behaviour, while the zonal meridional means show a more homogeneous behaviour. Here lower TCWV values are found in coincidence of longitudes where we have extended land presence in the equatorial region. The mean and the absolute standard deviation show similar features, with higher values in the ITCZ region. As can be noticed, the seasonal movements of the ITCZ from the north (in Northern Hemisphere summer) to south (in Northern Hemisphere winter) can be clearly detected (Castelli et al., 2018). In general, the standard deviations in the region where the TCWV maximum is located are of the order of 15 % and up to 20 % in polar regions. The zonal means reflect the shift of the ITCZ during the year, reaching maximum values of TCWV in Northern Hemisphere in July.

The quality of the AIRWAVEv2 dataset is evaluated through the same method adopted for AIRWAVEv1 and reported in Papandrea et al. (2018). The AIRWAVEv2 dataset is compared to the TCWV obtained from the SSM/I satellite and to data available from the ARSA.

In this contest, these two datasets are complementary, as the SSM/I TCWV are not retrieved for measurements in proximity to coasts (minimum distance of about 60 km), while the selected ARSA stations are located in coastal areas.

The SSM/I dataset has been produced from the SSM/I instrument series on board the DMSP polar orbit satellites since 1987. For this comparison, we used the 0.25^∘ v7 daily product obtained from the F13 satellite produced by remote-sensing systems. In fact, in contrast to the other DMSP satellites, the local time of the ascending node of F13 (18:00 UTC) is more stable than the one of the other satellites with only a variation of 1 h during the whole mission. For the comparison with SSM/I we used the AIRWAVEv2 data aggregated at 0.25^∘ resolution, covering the time period from 1995 to 2009. The ARSA dataset spans from January 1979 to the present and contains water vapour concentration profiles at specific pressure levels. To obtain the TCWV we vertically integrate these values. For the comparison with the ATSR data, only stations surrounded, even partially, by water are used. More details about the selection of ATSR and ARSA coincident data are reported in Papandrea et al. (2018).

Figure 5 Zonal means of TCWV for AIRWAVEv1 (blue), AIRWAVEv2 (red) and correlative measurements (black): SSM/I (a) or ARSA (b). The data have been averaged in 2^∘ latitude bins. AIRWAVE TCWV standard deviations are also reported.

The zonal means (calculated in bins of 2^∘ for all the datasets and reported in Fig. 5) show that AIRWAVEv2 data are of high quality with respect to both radiosondes and satellite data at all latitudes. In Fig. 5, the AIRWAVEv1 data are overplotted for comparison. The improvement in the performance of the new dataset is clearly visible at all latitudes, and in particular for regions at latitudes higher than 45–50^∘, where the negative values obtained with AIRWAVEv1 disappear. In addition, a significant reduction of the spread is highlighted.

Figure 6 AIRWAVE TCWV vs. SSM/I TCWV (a) or ARSA TCWV (b). The bin size is 2.5 kg m⁻². The colour scale indicates the number of elements of the histogram. The data cover the period from 1991 to 2012.

In Fig. 6 we show the bidimensional histograms of the comparisons between AIRWAVEv2 and SSM/I (panel a) and ARSA (panel b). The SSM/I measurements are homogeneously distributed over the globe, while ARSA radiosounding stations are mainly located at midlatitudes (see Fig. 7) and this is reflected in the bulk of the ARSA TCWV values ranging between 0 and 30 kg m⁻².

Figure 7 Average absolute (a) and relative (b) TCWV differences (SSM/I-AIRWAVE) at $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ spatial resolution. Average ARSA-AIRWAVE TCWV differences are overplotted with circles; the size of the symbols is proportional to the number of matches (see legend). The data cover the period from 1991 to 2012.

Globally, a good correlation is obtained from the datasets, as highlighted by the correlations (0.948 with SSM/I and 0.918 with ARSA) and bias values (0.02±4.79 kg m⁻² with respect to SSM/I and 0.19±6.12 kg m⁻² with respect to ARSA). We highlight that, in the validation exercise, we compare SSM/I data in coincidence with AIRWAVEv2 ones. Since AIRWAVE is applicable only to clear-sky measurements, this is a method to filter out SSM/I cloudy TCWV and thus to avoid biases due to different sensitivities related to the used spectral range. When comparing AIRWAVE data to radiosondes, the small bias we found demonstrates that AIRWAVE TCWV are also sensitive to low atmospheric levels. In particular, the comparison with the same histograms of AIRWAVEv1 (see Fig. 1 of Papandrea et al., 2018) highlights the correction of negative values in polar and coastal regions in the new version of the dataset.

Table 2 AIRWAVE TCWVs compared with SSM/I and ARSA stations. Results are also given for AIRWAVEv1. Absolute differences along with the standard deviations are reported for the global (all latitudes), equatorial (25^∘ S–25^∘ N), midlatitude (25–60^∘ S and 25–60^∘ N) and polar (> 60^∘ N or > 60^∘ S) scenarios. Average values for ATSR-1, ATSR-2 and AATSR are also provided.

Figure 7 reports the geographical distribution of the mean TCWV differences with respect to SSM/I and ARSA in absolute and percentage values. In comparison to AIRWAVEv1 (see Fig. 3 of Papandrea et al. (2018)) the differences with SSM/I are reduced at all latitudes. The longitudinal patterns of the differences are similar to the ones of AIRWAVEv1, except for the equatorial pacific region where AIRWAVEv2 shows a slightly higher positive bias. The reasons for this behaviour are under investigation. In the majority of coastal regions, where AIRWAVEv1 underestimated the TCWV, AIRWAVEv2 is now in agreement with ARSA results (no SSM/I data close to the coast).

In Table 2 we summarise the results of this comparison for both ARSA and SSM/I and for different scenarios and missions. The average bias is about 0.0±4.7 kg m⁻² with respect to SSM/I and 0.2±6.1 kg m⁻² with respect to ARSA. If we compare these results with the ones for AIRWAVEv1 (reported in Table 2 to ease the comparison) we can clearly see the improvement in both the biases and the standard deviations (0.7±5.7 kg m⁻² with respect to SSM/I and 0.8±7.7 kg m⁻² with respect to ARSA). As can be seen from Figs. 5 and 7 and from the results in Table 2 the improvement in the bias is obtained at all latitudes; it is, however, more evident in polar regions (from 5.5±5.1 kg m⁻² in AIRWAVEv1 to 1.3±3.5 kg m⁻² in AIRWAVEv2 versus SSM/I and from 4.1±6.5 kg m⁻² in AIRWAVEv1 to 0.9±4.6 kg m⁻² in AIRWAVEv2 when using ARSA). Slight differences between the three ATSR missions are consistent with the related uncertainties.

Figure 8 SSM/I-AIRWAVE TCWV (orange) and ARSA-AIRWAVE TCWV (red) monthly mean trends. The difference between the correlative measurements and AIRWAVE TCWV ± SD is also reported.

In Fig. 8 we show the monthly mean evolution of the differences (and their standard deviation) between the TCWV obtained from correlative measurements and AIRWAVEv2. As for AIRWAVEv1, these differences are quite stable over time, with the exception of the beginning of the ATSR-1 mission (1991–1994). As explained in Papandrea et al. (2018), this could be due to the failure of a 3.7 µm channel that impacted the ATSR cloud screening. In general, the differences with respect to the radiosonde exhibit a higher seasonality due to pronounced variability of atmospheric and surface conditions in coastal areas.

It is worth noticing that the spread of AIRWAVEv2 is always smaller than the one of AIRWAVEv1. This is partially due to the new algorithm, which reduces the random error component due to noise in the retrieved TCWV.

The above-described results indicate that the AIRWAVEv2 algorithm reduces the global bias with respect to SSM/I, from about 0.7 kg m⁻² of AIRWAVEv1 to 0.0 kg m⁻² and improves the standard deviation (SD) of up to 20 % with respect to AIRWAVEv1. In AIRWAVEv2, the SD values are essentially constant for all the scenarios and missions, highlighting the unbiased nature of the dataset with respect to SSM/I (cloud-free). When using ARSA data the bias reduces from 0.8 to 0.2 kg m⁻² and the standard deviations by 21 %.

4 Discussion and conclusions

The second version of the AIRWAVE TCWV dataset described in this work has been validated against ARSA and SSM/I equivalent products.

As expected also from the analysis of synthetic retrievals, the most significant AIRWAVEv2 improvement is achieved at polar latitudes. In polar regions the bias versus SSM/I improves of 4.2 and of 3.2 kg m⁻² versus ARSA. In both cases the standard deviations are reduced of about 1.6–1.9 kg m⁻². However, improvements at midlatitudes are also found. The average bias with respect to SSM/I improves by about 0.7 kg m⁻² and the standard deviation is reduced by about 1 kg m⁻². In the case of validation against radiosondes, the bias in AIRWAVEv2 is reduced by about 0.6 kg m⁻² with respect to AIRWAVEv1 and the standard deviation is reduced by 1.6 kg m⁻².

These improvements are due to two factors. First the AIRWAVEv2 retrieval parameters now account for the atmospheric variability. Secondly the implementation of the new algorithm explicitly takes into account the geometry and latitude dependence of each pixel, allowing possible artefacts due to approximations and a posteriori corrections to be overcome. No statistically significant trend can be found in the comparison with SSM/I and ARSA in both versions of the database, while a seasonal dependence of the differences is observed, with a larger bias in July and August, mainly due to the differences in midlatitude north TCWV retrievals. In general, we find slightly drier results with respect to ARSA and SSM/I with both versions. This is possibly due to the fact that the temporal mismatch between the ATSR and the correlative measurements does not allow all SSM/I TCWV retrieval obtained under cloudy conditions to be excluded or cloud masks to be wrongly assigned. As discussed, the use of retrieval parameters that are calculated in conditions different from the ones present in the observed scenario can cause biases in the obtained TCWV. We point out that the major source of errors in the retrieved TCWV comes from the temperature profile assumptions, while erroneous assumptions of other gases (e.g. HNO₃, CFC-11, CFC-12, CO₂) have an almost negligible impact. The obtained RMSE value of about 7 % is of the same order of this error. The error quantification given in this work allows the users to get a better insight of the AIRWAVEv2 TCWV dataset and its related quality.

Besides the improvements on ATSR TCWV retrievals given by the AIRWAVEv2 dataset, the method described in this work can be the basis for a similar approach for SLSTR (Sea and Land Surface Temperature Radiometer, on board Copernicus Sentinel-3, Donlon et al., 2012).

Appendix A: Evaluation of the main systematic and random errors in AIRWAVE version 2

This appendix provides an estimate of the main error sources (random and systematic) that affect the AIRWAVEv2 dataset. To facilitate the use of the errors together with the TCWV contained in the dataset, we summarise the results of this Appendix in Table A1. AIRWAVEv2 retrieval parameters account for atmospheric and surface variability of the observed scenes. However, real observed scenes can obviously deviate from the simulated ones. An evaluation of the errors induced by these deviations is then required. For this reason we analyse, through the use of synthetic radiances, the major sources of errors that can affect the AIRWAVEv2 retrieval, i.e. (in order of importance) the influence of the atmospheric temperature profile variations, the impact of sea emissivity changes due to the wind and the impact of interfering atmospheric species. Furthermore, an estimate of the random error component due to the noise is also reported.

A1 Retrieval approximations

In order to evaluate the impact of systematic errors due to the retrieval approximations, we ran some tests of TCWV retrievals from simulated radiances. We used top-of-atmosphere (TOA) subsatellite track brightness temperatures (BTs) simulated with the RTM at 10.8 and 12 µm in nadir and forward geometries as input to the retrieval chain (no random noise was added). The BTs were produced using different TCWV amounts for each given scenario (e.g. water vapour profiles were multiplied by 0.5, 0.75, 1.25 and 1.5), while all the other atmospheric profiles were kept constant. The results of this exercise show that TCWV is correctly retrieved with an error of ±3 %. For comparison purposes, we performed a similar test using the AIRWAVEv1 approach. For AIRWAVEv2 tests, the same atmospheric conditions used to compute the retrieval parameters have been adopted, while for AIRWAVEv1 larger errors are expected due to the fact that the atmospheric variability is not taken into account. Despite this, this analysis shows that AIRWAVEv1 performs well for medium-high TCWV (> 30 kg m⁻²), where the differences are below 10 %. For TCWV values between 10–20 kg m⁻², the differences range between −10 % and −30 %. AIRWAVEv1 underestimates the low TCWV values (below 10 kg m⁻²) with differences from −50 % to $-\mathrm{150}/-\mathrm{200}$ %. This was also reported in Casadio et al. (2016), where the authors compared AIRWAVE TCWV with collocated ECMWF counterparts, showing a dry bias at high latitudes (where the TCWV values are small). The tests on synthetic radiances indicate that the AIRWAVEv2 parameterisation solves this issue.

Figure A1 CO₂, HNO₃, CFC-11 and CFC-12 spectral lines and ATSR filter functions (arbitrary units).

A2 Atmospheric temperature and water vapour profiles and SSTs

One of the main error sources that can affect the AIRWAVE TCWV retrieval is the assumption of a fixed temperature profile. Actually, the atmospheric opacity is closely linked to the atmospheric density and thus to the atmospheric temperature. To estimate the impact of temperature on the retrieved TCWV, 20 different temperature profiles were randomly perturbed by ±3 K on a 1 km equispaced altitude grid. In order to also account for possible changes related to SST variation we change the SST accordingly to the value of the temperature in the lowest layer. Then, simulated BTs were produced with the RTM and were used to perform the TCWV retrievals for the three instruments in equatorial, midlatitude and polar July conditions for the Northern Hemisphere and the results were compared with respect to the unperturbed case. In the third column of Table A1 we summarise the findings of this analysis for each of the three instruments, reporting the SD of the difference both in absolute and in percentage values. The impact of these perturbations is of the order of 6 % and is higher in the equatorial and midlatitude regions and lower at the poles. Indeed, in the equatorial region, due to the higher water vapour content, the atmosphere is more opaque than at the poles so that temperature variations have a larger impact on the retrieved TCWV. These tests were also run by varying the atmospheric profile alone. Very similar but slightly higher errors are found in these cases. Another relevant error sources can be due to differences in the water vapour profile shape. To evaluate this error, we varied the water vapour profile randomly up to 5 % at each atmospheric level. At maximum the impact is 1 % in the equatorial case (atmospheric opacity; see the fourth column in Table A1).

Table A1 Impact of different error sources on AIRWAVEv2 dataset.

A3 Wind speed

A further source of error that can affect the TCWV retrievals is the value used for the sea emissivity, which directly enter Eqs. (10) and (11). Sea emissivity depends on wavelength, sea surface temperature, viewing angles and wind speed. As stated in Sect. 2.2, in AIRWAVEv2 we accounted for emissivity variations due to the viewing angles and sea surface temperatures. All the calculations were made at a fixed wind value of 3 m s⁻¹. In order to assess the possible systematic effects due to wind variations on the retrieved values, we varied the emissivity according to the wind speed at three values as tabulated in the University of Edinburgh emissivity database: 1, 10 and 25 m s⁻¹.

The emissivity has a different value and spectral behaviour with different wind speeds. The fourth column of Table A1 reports the difference in the TCWV retrieved for simulated measurements with the 25 m s⁻¹ wind speed with respect to the reference case. As expected, the impact is almost negligible in the equatorial band, where the higher opacity of the atmosphere reduces the sensitivity of ATSR measurements to the surface conditions, while increasing it toward the poles where, due to the low atmospheric opacity, the surface effects become relevant with respect to the atmospheric component. Furthermore, in the case of polar conditions, the effects of the wind on the retrieved TCWV are not linear, with enhanced variations for wind speed of 25 m s⁻¹, as not only the intensity but also the spectral shape of the surface emissivity varies in function of the wind speed. The possible presence of white caps, generated by high speed winds, has not been considered in this study.

A4 Interfering atmospheric constituents

The AIRWAVE algorithm, in both versions, accounts for the contribution to the radiance of the two main gases active in the ATSR channels (H₂O and CO₂). However some other species have spectroscopic features in the ATSR channels. In Fig. A1, the CO₂, HNO₃ and CFCs spectra in the 10–13 µm wavelength range are shown, along with the ATSR filter functions (all in arbitrary units). In order to have a complete view of the possible error components, we assessed the impact of interfering species in case their abundance differs from the one used in the reference scenarios to compute the retrieval parameters.

Among the considered species, HNO₃ shows significant latitudinal and seasonal variability, while CO₂ and CFCs exhibit interannual trends (Remedios et al., 2007). For this reason we separately accounted for the effects due to latitudinal and interannual variability. We used the IG2 database version 4.1 and our RTM to produce synthetic BTs. For each season and latitude band we generated synthetic BTs using the proper IG2 atmospheric status, except for the profile of the investigated species, for which we used all the different available profiles. The generated BTs were then used to retrieve the TCWV to assess the systematic error induced by the expected variability of the interfering species.

Latitudinal and seasonal variations of HNO₃ impact ATSR-1 and ATSR-2 BTs more than AATSR, because of the different shapes of the TIR filter functions, and can produce differences up to 0.3 K in the 11 µm band (midlatitude vs. tropical north in January). For CFC-11 seasonal differences in the tropics are of the order of 0.001 K in the 12 µm band and for CFC-12 we get 0.03 K in the 11 µm band, while latitudinal variations reach 0.04 K from tropical to midlatitude atmospheres for both channels. CO₂ latitudinal variations can produce a maximum difference of 0.003 K on nadir BT. The impact of maximum latitudinal variations of HNO₃, CFCs and CO₂ on the retrieved TCWV is reported in Table A1. The largest contribution is due to HNO₃ latitudinal variation and is of the order of 0.6 %, while the CFCs latitudinal variations produce differences of 0.15 % and CO₂ of 0.01 % only. Furthermore, the impact of using mid latitude profiles instead of tropical profiles for all species but H₂O has a maximum impact of 0.22 %. Seasonal variations are almost negligible for CO₂ and CFCs, while they are only 0.6 % for HNO₃. Thus we can safely assume that the latitudinal and seasonal variations of the interfering species represent a minor error source for the AIRWAVE TCWV for both AIRWAVEv1 and AIRWAVEv2.

To also evaluate the impact of interannual variations for all the ATSR series, we would need the CO₂ and CFCs profiles from 1991 to 2012. However, as mentioned in the previous sections, the IG2 database contains data from 2002 onwards. In this work, the CFCs and CO₂ profiles from 1991 to 2001 were inferred by scaling the 2002 profiles using the trend given in the last IPCC report for CO₂ and on the Mauna Loa observatory website (Global Monitoring Division, 2018 website; see also Aoki et al., 2003 and Minschwaner et al., 2013) for CFC-11 and CFC-12. The impact on the ATSR BTs due to the interannual variations of CO₂ and CFCs has been evaluated for each mission. We calculated the synthetic spectra using the CFC-11, CFC-12 or CO₂ profile for the initial and final year of each mission. Then we calculated, for each instrument, the differences in retrieved TCWV at the beginning and at the end of each mission. The results for CO₂ and CFCs are shown in Table A1. The influence of CO₂ annual variations on the retrieved TCWV is of the order of 0.004 kg m⁻² (< 0.01 %). For CFCs we obtain 0.002 kg m⁻². We can then conclude that the impact of VMR latitudinal and seasonal variations on retrieved TCWV is negligible (maximum value 0.6 %), and that the systematic effect of CO₂ and long-term CFCs variations over the missions are even smaller (0.07 %).

A5 Noise

Finally we analyse the impact of the measurement noise on the retrieved TCWV (see last column of Table A1). The measurement noise was simulated by applying a random perturbation of ±0.037 K to the BTs of the two channels of ATSR-2 and AATSR and a perturbation of ±0.1 K on the ATSR-1 channels (Smith et al., 2012). For each instrument, we have generated 1000 values and we have evaluated the standard deviation of the obtained TCWV that we report in absolute and percentage values in Table A1. The standard deviation is maximum at the poles, as expected, since there the TCWV and S ∕ N ratio are lower. For AIRWAVEv2 we get 18 % for ATSR-1 and 6 % for ATSR-2 and AATSR in the worst case. Note that AIRWAVEv2 approach has, for all scenarios, better performance with respect to AIRWAVEv1 (see Casadio et al., 2016).