A cloud algorithm based on the O₂-O₂ 477 nm absorption band featuring an advanced spectral fitting method and the use of surface geometry-dependent Lambertian-equivalent reflectivity

We discuss a new cloud algorithm that retrieves an effective cloud pressure, also known as cloud optical centroid pressure (OCP), from oxygen dimer (O₂-O₂) absorption at 477 nm after determining an effective cloud fraction (ECF) at 466 nm, a wavelength not significantly affected by trace-gas absorption and rotational Raman scattering. The retrieved cloud products are intended for use as inputs to the operational nitrogen dioxide (NO₂) retrieval algorithm for the Ozone Monitoring Instrument (OMI) flying on the Aura satellite. The cloud algorithm uses temperature-dependent O₂-O₂ cross sections and incorporates flexible spectral fitting techniques that account for specifics of the surface reflectivity. The fitting procedure derives O₂-O₂ slant column densities (SCDs) from radiances after O₃, NO₂, and H₂O absorption features have been removed based on estimates of the amounts of these species from independent OMI algorithms. The cloud algorithm is based on the frequently used mixed Lambertian-equivalent reflectivity (MLER) concept. A geometry-dependent Lambertian-equivalent reflectivity (GLER), which is a proxy of surface bidirectional reflectance, is used for the ground reflectivity in our implementation of the MLER approach. The OCP is derived from a match of the measured O₂-O₂ SCD to that calculated with the MLER method. Temperature profiles needed for computation of vertical column densities are taken from the Global Modeling Initiative (GMI) model. We investigate the effect of using GLER instead of climatological LER on the retrieved ECF and OCP. For evaluation purposes, the retrieved ECFs and OCPs are compared with those from the operational OMI cloud product, which is also based on the same O₂-O₂ absorption band. Impacts of the application of the newly developed cloud algorithm to the OMI NO₂ retrieval are discussed.

1 Introduction

Satellite ultraviolet and visible (UV–vis) nadir backscattered sunlight trace-gas algorithms need accurate estimates of cloud parameters in order to produce high-quality data sets. Because of complexity of cloud effects on the radiation field in the atmosphere, clouds in trace-gas algorithms are treated using multiple simplifying assumptions. Among them, the fundamental assumptions are (1) the independent pixel approximation (IPA) that neglects the horizontal transport of radiative energy between the clear-sky and overcast subpixels, and (2) the assumption of horizontally and vertically homogeneous clouds that substantially simplifies radiative transfer in the clouds. For trace-gas retrievals it is important to estimate the photon path lengths in the atmosphere that determine trace-gas absorption and thus affect the measured top-of-atmosphere (TOA) radiance. The photon path lengths in a cloudy atmosphere are determined by the following most important cloud parameters, the geometrical cloud fraction and the cloud vertical extinction profile (Stammes et al., 2008). Other cloud parameters: the cloud phase, the cloud particle shape, and the particle size distribution, that determine the cloud phase scattering function are usually not considered. Because of the limited informational content of TOA radiances, all of those parameters cannot be independently retrieved from the radiance measurements. That is why the additional cloud assumptions should be adopted. For instance, if we retrieve the cloud parameters from an oxygen absorption band and assume a model of scattering cloud (Loyola et al., 2018), we have to adopt a priori values of the cloud microphysical parameters and cloud vertical extent assuming a homogeneous cloud layer and to add information to cloud fraction from other measurements.

One of the simplest cloud models is the so-called mixed Lambertian-equivalent reflectivity (MLER) model. Many trace-gas algorithms are based on the MLER model. For example, the MLER model is currently used in most trace-gas (Bucsela et al., 2013; Veefkind et al., 2006) and cloud (Acarreta et al., 2004; Joiner and Vasilkov, 2006; Veefkind et al., 2016) retrieval algorithms for the Ozone Monitoring Instrument (OMI) (Levelt et al., 2006), a Dutch–Finnish UV–vis sensor flying on the National Aeronautics and Space Administration (NASA) Aura satellite. For each field of view (FOV) the MLER model treats cloud and ground as horizontally homogeneous, opaque Lambertian surfaces and mixes them using IPA. According to the IPA, the measured TOA radiance is a sum of the clear-sky and overcast subpixel radiances that are weighted with an effective cloud fraction (f) (Stammes et al., 2008). f is calculated by inverting the equation

$$\begin{array}{}\text{(1)}& I_{\mathrm{m}}=I_{\mathrm{g}}\left(R_{\mathrm{g}}\right)(\mathrm{1}-f)+I_{\mathrm{c}}\left(R_{\mathrm{c}}\right)f\end{array}$$

at a wavelength not substantially affected by rotational Raman scattering (RRS) or atmospheric absorption, where I_m is the measured TOA radiance, I_g and I_c are the precomputed clear-sky (ground) and overcast (cloudy) subpixel TOA radiances, and R_g and R_c are the corresponding ground and cloud Lambertian-equivalent reflectivities (LERs), respectively.

The MLER model typically assumes R_c=0.8 (Koelemeijer et al., 2001; McPeters et al., 1996). This assumption more accurately accounts for Rayleigh scattering in partially cloudy scenes (Ahmad et al., 2004; Stammes et al., 2008) and also accounts for scattering/absorption that occurs below a thin cloud. In this paper we also adopt R_c=0.8 for the O₂-O₂ cloud algorithm. The ground reflectivity is usually taken from an LER or surface albedo climatology derived from satellite observations (Herman and Celarier, 1997; Kleipool et al., 2008). In reality, reflection of incoming direct and diffuse solar light from land or ocean surfaces is sensitive to the sun–sensor geometry. This dependence is described by the bidirectional reflectance distribution function (BRDF). Impact of surface BRDF on the cloud and trace-gas retrievals has been studied since 2010. For instance, Zhou et al. (2010) reported that accounting for surface BRDF effects can change NO₂ retrievals by up to 20 %. The newest advances in this field and more references can be found in the latest paper by Lorente et al. (2018). To account for the BRDF, we developed a new model of geometry-dependent Lambertian equivalent reflectivity (GLER) that was implemented within the existing OMI cloud and NO₂ retrieval algorithms (Vasilkov et al., 2017). This implementation only required changes to the input surface reflectivity database, thus simplifying the introduction of BRDF effects.

The MLER model compensates for photon transport within a cloud by placing the Lambertian surface somewhere in the middle of the cloud instead of at the top (Vasilkov et al., 2008). As clouds are vertically inhomogeneous, the pressure of this surface does not necessarily correspond to the geometrical center of the cloud, but rather to the so-called optical centroid pressure (OCP) (Joiner et al., 2012; Sneep et al., 2008; Vasilkov et al., 2008). The cloud OCP can be thought of and modeled as a reflectance-averaged pressure level reached by backscattered photons (Joiner et al., 2012). Cloud OCPs are the appropriate quantity for use in trace-gas retrievals from satellite instruments (Joiner et al., 2006, 2009; Vasilkov et al., 2004). Cloud-top pressures derived from thermal infrared (IR) measurements are not equivalent to OCPs and do not provide good estimates of solar photon path lengths through clouds that are needed for trace-gas retrievals from UV–vis backscatter measurements (Joiner and Vasilkov, 2006; Joiner et al., 2012; Vasilkov et al., 2008).

The main goal of this paper is to document new approaches to the development of a cloud algorithm based on the O₂-O₂ absorption band at 477 nm (Yang et al., 2015). These approaches include an advanced spectral fitting algorithm for TOA radiances and the use of surface GLERs to replace climatological LER data sets. This new cloud algorithm is intended for the operational OMI NO₂ algorithm and is planned as a backup cloud algorithm for the upcoming Tropospheric Emissions: Monitoring Pollution (TEMPO) geostationary mission (Zoogman et al., 2017). Our spectral fitting procedure is similar to that developed by Marchenko et al. (2015) for NO₂ slant column density (SCD) retrieval. It relies on the temperature-dependent O₂-O₂ cross sections (Thalman and Volkamer, 2013) and derives the O₂-O₂ SCD using O₃, NO₂, and H₂O slant column estimates from independent OMI algorithms. This is an implementation choice that is designed to minimize potential errors due to interference among O₃, NO₂, and O₂-O₂ spectral absorption features during the fitting procedure. We apply the new cloud algorithm to the OMI NO₂ retrievals and show NO₂ column changes related to the use of the new cloud algorithm.

2 Data and methods

2.1 OMI and MODIS data

We use several data sets from the OMI and Moderate resolution Imaging Spectroradiometer (MODIS) instruments flying on the NASA Aqua and Terra satellites. OMI is a spectrometer that acquires Earth and solar spectra at UV–vis wavelengths from 270 to 500 nm with a spectral resolution of approximately 0.5 nm. The OMI ground footprint varies; near the nadir, it is approximately 12 km along the satellite track and 24 km across the 2600 km track. The footprint size increases towards the swath edge. We use TOA radiance and solar irradiance in the OMI vis channel to retrieve cloud parameters and NO₂ amounts.

The MODIS-derived BRDF kernel coefficients from the 16-day MCD43GF data set (Schaaf et al., 2011) are used to compute GLERs over land for the OMI swath (Vasilkov et al., 2017). The kernel coefficients are provided for snow-free land and permanent ice at a high spatial resolution. Over transient snow-covered regions, we retain the standard climatological LER of Kleipool et al. (2008) that was routinely used for the previous cloud retrievals.

2.1.1 GLER computation

The BRDF kernel coefficients are averaged over an OMI FOV and used to calculate the TOA radiance for a given observational geometry assuming pure Rayleigh scattering in the atmosphere. For radiative transfer (RT) calculations, we use the vector linearized discrete ordinate radiative transfer (VLIDORT) code (Spurr, 2006). VLIDORT computes the Stokes vector in a plane-parallel atmosphere with a Lambertian or non-Lambertian underlying surface. It has the ability to deal with attenuation of solar and line-of-sight paths in a spherical atmosphere, which is important for large solar zenith angles (SZAs) and viewing zenith angles (VZAs). VLIDORT accounts for polarization at the ocean surface using a full Fresnel reflection matrix.

The TOA radiance computed by VLIDORT is then inverted to derive GLER using the following exact equation:

$$\begin{array}{}\text{(2)}& I_{\mathrm{TOA}}=I_{\mathrm{0}}+{\displaystyle \frac{RT}{\mathrm{1}-R{S}_b}},\end{array}$$

where I₀ is the TOA radiance calculated for a black surface, R is the GLER, T is the total (direct+diffuse) solar irradiance reaching the surface converted to the ideal Lambertian-reflected radiance (by dividing by π) and then multiplied by the transmittance of the reflected radiation between the surface and TOA in the direction of a satellite instrument, and S_b is the diffuse flux reflectivity of the atmosphere for the case of its isotropic illumination from below (Dave, 1978). All quantities, I₀, T, and S_b are calculated using a known surface pressure. We use a monthly climatology of surface pressure taken from the Global Modeling Initiative (GMI) chemistry transport model driven by the NASA Global Modeling and Assimilation Office (GMAO) GEOS-5 global data assimilation system (Rienecker et al., 2011) with a spatial resolution of 1^∘ latitude by 1.25^∘ longitude. Surface pressure for each OMI pixel P_s is calculated as follows:

$$\begin{array}{}\text{(3)}& P_{\mathrm{s}}=P_{\mathrm{s}}\left(\mathrm{GMI}\right)\exp (-\mathrm{\Delta }z/H),\end{array}$$

where P_s(GMI) is the GMI surface pressure at a resolution of $\mathrm{1}{}^{\circ }\times \mathrm{1.25}{}^{\circ }$, $\mathrm{\Delta }z=z-z\left(\mathrm{GMI}\right)$, z is the terrain height of the OMI pixel from a digital elevation model, z(GMI) is the terrain height at resolution of $\mathrm{1}{}^{\circ }\times \mathrm{1.25}{}^{\circ }$, $H=\left(kT\right)/\left(Mg\right)$ is the scale height, where k is the Boltzmann constant, T is the GMI air temperature at the surface, M is the mean molecular weight of air, and g is the acceleration due to gravity.

To calculate TOA radiance over water surfaces, we account for both light specularly reflected from a rough water surface and also for diffuse light backscattered by water bulk and transmitted through the water surface. Reflection from the water surface is described by the Cox–Munk slope distribution function as implemented in Mishchenko and Travis (1997). Diffuse light from the ocean is calculated using a Case 1 water model that has chlorophyll concentration as a single input parameter. Bidirectionality of the underwater diffuse light is accounted for following Morel and Gentili (1996).

More details about the GLER computation can be found in Vasilkov et al. (2017). An important update of our ocean surface model is the use of a variable wind speed instead of a single climatological wind speed of 5 m s⁻¹ as in Vasilkov et al. (2017). Retrievals of wind speed are taken from the Advanced Microwave Scanning Radiometer for the Earth Observing System (AMSR-E) that flies on NASA's Aqua satellite (with Aura-OMI closely following Aqua). The use of wind speed from the AMSR-E measurements improves in the GLER over ocean. Thanks to the higher spatial resolution of AMSR-E, it is possible to match fine structure of the wind field to the TOA radiances and GLERs over the sun-glint-affected areas.

2.2 The O₂-O₂ slant column density fitting algorithm

The operational OMI O₂-O₂ SCD retrieval (Acarreta et al., 2004; Veefkind et al., 2016) uses the differential optical absorption spectroscopy (DOAS) (Platt and Stutz, 2006) approach, simultaneously retrieving SCDs of O₂-O₂ and O₃ in the 460–490 nm wavelength interval and using single-temperature O₂-O₂ cross sections and a first-degree polynomial approximating the wavelength dependence of reflectances in the fitting window.

Here, we generally follow the approach developed by Marchenko et al. (2015) for the NO₂ SCD estimates. Instead of simultaneous retrieval of coefficients of multiple parameters as takes place in the classical DOAS formalism, we divide the problem into a series of sequential steps (Fig. 1).

Figure 1 Flow diagram of the O₂-O₂ SCD retrieval algorithm. The algorithm input comprises the OMI monthly mean solar irradiances; the radiances (dependent on wavelength, line of sight (row), and position (along orbit)); the laboratory cross sections of O₃, NO₂, and H₂O (X sections); the atmospheric (RS air) and liquid-water (RS water) Raman scattering spectra (all X sections convolved with the row- and wavelength-dependent OMI instrument line-shape functions); and the OMI cloud-fraction (CF) estimates provided by an independent retrieval. RS denotes the amplitudes of the combined air and water Raman scattering spectrum.

Step 1 involves the removal of interfering trace-gas absorption. The spectral range chosen for the O₂-O₂ SCD retrievals is affected by relatively strong O₃ absorption that, in most cases, distorts the O₂-O₂ profiles (Fig. 2). The same applies to NO₂ absorption over polluted regions (e.g., the Beijing area; see Fig. 3) and, to a far lesser extent, the mainly equatorial regions over the open-water Pacific, where the H₂O absorption may distort the flanks of the broad O₂-O₂ profiles (Fig. 2). Note the clear presence of the ozone feature around λ∼462 nm, as well as the large distortion of the O₂-O₂ profile caused by the broad ozone absorption around λ∼482 nm. In this particular example, the only easily recognizable Ring spectrum feature is seen at λ∼486.5 nm. The gradual ∼13 % change in reflectances between 450 and 500 nm comes from a combination of comparable strength signals: the ozone absorption and the Rayleigh-scattering component. For better guidance, in Figs. 2–5 we show scaled absorption spectra of the main trace gases that may contribute to the general appearance of a reflectance spectrum. We keep the same plotting style for Figs. 2–5, though noting that in each particular example one may see quite a different impact from the same absorption constituent. For example, while the NO₂-related signal barely registers in Fig. 2 (practically unpolluted region), the heavily contaminated Beijing area (Fig. 3) shows a clear presence of the NO₂ absorption at λ∼457–466 nm (three features), λ∼475 and 480 nm (these two are superposed on the broad O₂-O₂ absorption), and the well-defined NO₂ absorption at λ∼489 nm that rivals the strength of the retrieved O₂-O₂ feature. In this particular case of the heavily NO₂-polluted region, the omnipresent O₃ absorption plays a far less important role compared to the spectrum shown in Fig. 2. In addition, though to a far lesser extent (when compared to ozone, Ring, and NO₂), the H₂O absorption may distort the flanks of the broad O₂-O₂ profiles, mainly in the equatorial regions over the open-water Pacific.

Figure 2 (a) Reflectance normalized at λ=464 nm (bold black line) for the OMI orbit 7921 from 10 January 2006, with row 14 and orbital exposure 1550 at 65.37^∘ N and 88.58^∘ E (high slant column ozone values). For reference, the arbitrarily shifted and scaled absorption spectra of H₂O (thin blue line), NO₂ (green), O₂-O₂ (red), and O₃ (cyan) are plotted in the upper portion of the panel. The arbitrarily scaled and shifted Ring patterns (as seen in reflectances) are shown in black. (b) The rectified O₂-O₂ absorption profile (i.e., the ratio of the data denoted by the red and blue lines in a, with additional adjustments to the blue-line data – see text for more details). The dashed black line shows the 1.0 reference level. In panel (b), the vertical dotted lines denote the wavelength range used in the SCD fits of the O₂-O₂ absorption profile.

The spectral domain chosen for the O₂-O₂ retrieval is not optimal for simultaneous O₃, NO₂, or H₂O estimates. Optimal fitting windows are the ∼290–340, 400–465, and 435–450 nm intervals for O₃, NO₂, and H₂O, respectively. Hence, to minimize the biases that may be introduced by the sub-optimal DOAS estimates of the interfering trace-gas species, we use the SCDs provided by independent OMI products: NO₂ and H₂O from OMNO2SCD (Marchenko et al., 2015) and O₃ from OMDOAO3 (Veefkind et al., 2006), and we remove the corresponding absorption features from the observed radiances. We find that, as expected, at large SZAs the corrections based on the UV O₃ SCD retrievals result in large spectral residuals pointing to systematic underestimates of O₃ absorption strength. This stems from the notion that the relatively (to the visual range) higher Rayleigh optical depth effectively masks the lower-atmosphere O₃ absorption. Hence, at $\mathrm{SZA}>\mathrm{80}{}^{\circ }$ we adjust the UV SCDs by a constant 1.25 coefficient. This helps to reduce the spectral residuals related to the underestimated O₃ absorption to a manageable (on average <0.1 %) level.

Step 2 closely follows the approach from Marchenko et al. (2015), comprising the simultaneous, iterative wavelength adjustment and Ring spectrum removal. At each FOV (row, 60 in total) the reflectances are produced from the individual, prefiltered earthshine radiances (from Step 1) normalized by the monthly-averaged OMI irradiances. These irradiances are iteratively adjusted to accommodate slight relative (radiances vs. irradiances) wavelength shifts. For estimates of the line-filling factors (i.e., the Raman scattering amplitudes), we use an appropriate combination of the air and water Raman scattering spectra (Vasilkov et al., 2002). We split the retrieval region into two “micro-windows”, 451–469 and 483–496 nm, and iteratively evaluate the wavelength shifts and the Raman-spectrum amplitudes in each window. For the final removal of the Raman scattering patterns, we use an average of the two micro-window estimates. The individual micro-window wavelength shifts are used for wavelength adjustments of irradiances in each micro-window, interpolating these estimates in the 469–483 nm domain occupied by the main O₂-O₂ absorption.

Figure 3 Similar to Fig. 2 but for row 44 and orbital scan number iTime=1315 pointing to the Beijing area. In (a), the thick blue line follows the reflectances after removal of the trace-gas (O₃ and H₂O) absorption, however with the Ring spectrum features remaining intact. This is to be compared to the adjacent black line that follows the original reflectances. Panel (b) shows the normalized O₂-O₂ profile as used in the SCD retrieval.

Figure 4 Similar to Fig. 3 but for the open-water nearly cloud-free (f<0.05) region of the Indian Ocean (54.02^∘ S, 106.91^∘ E, OMI orbit 7791, 1 January 2006). The thick blue line follows these reflectances after removal of the Ring patterns and the trace-gas (O₃, NO₂, and H₂O) absorption. The red line shows the piecewise fit to the blue line.

Figure 5 Similar to Fig. 4 but for the Sahara desert (OMI orbit 8013, 16 January 2006, row 20, orbital exposure 1180). Panel (b) shows the rectified O₂-O₂ absorption profile, i.e., the ratio of the data denoted by the red and blue lines in (a), with additional adjustments to the blue-line data – see text for more details.

Step 3 involves normalizing the O₂-O₂ profile in preparation for SCD evaluations. We deem this step to be the most important procedure; it may change the outcome by as much as ∼20 % in extreme cases such as the open-water scenes (Fig. 4) and Sahara desert (Fig. 5) representing two extremes and the remaining cases falling in between. Figures 4a and 5a show the observed reflectances before (black lines) and after (blue lines) the removal of trace-gas absorption and the Raman line-filling patterns, and the red lines follow the adopted continuum fits. Note the profound difference among the wavelength dependencies of the reflectances in these extreme cases. While the cloud-free, open-water case (Fig. 4) is predominantly Rayleigh controlled, leading to a steep decline in reflectances, the much brighter Sahara surface controls the appearance of the radiances at long OMI wavelengths, leading to the gradual increase in reflectances (Fig. 5). Figures 4b and 5b show normalized O₂-O₂ profiles. In a case-by-case study of the presumably cloud-free areas, we have found that various combinations of linear functions fitting the flanks of the O₂-O₂ profile lead to gross underestimates (mainly over open-water areas) or overestimates (deserts and semideserts) of the retrieved scene pressures that are directly linked to the biases in the SCD evaluations. Hence, we have implemented a more flexible approach, defining two broad categories of the surface reflectances and applying different fitting approaches to each of them.

The reflectances from Step 2 are averaged in 2 nm intervals, providing a set of estimates at λ=463 and 495 nm that are partitioned into two general categories. The first broad category comprises all the relatively cloud-free low-reflectance scenes, with $\frac{r\left(\mathrm{463}\right)}{r\left(\mathrm{495}\right)}>\mathrm{1.05}$ and r(463)<0.25. The second class includes the remaining scenes. For both categories, the fitting starts from applying the third-degree polynomial to the 459–466 and 484–494 nm regions, identifying and eliminating large (∼5σ, i.e., ±0.5 %) deviations and then repeating the procedure, ultimately normalizing the reflectances in the 450–500 nm range by the fit. The normalized reflectances in the 459–465 and 484–490 nm intervals are refitted with a first-degree polynomial and, again, all reflectances in the 450–500 nm range are renormalized by this fit. This concludes the fitting for the second category of scenes. However, at this point the fitting proceeds for the first class of the relatively cloud-free, low-reflectance scenes. Yet again, the normalized radiances in the 465–470 and 482–487 nm intervals are iteratively (rejecting the large ±1 % deviations) refitted with a second-degree polynomial; then this fit is applied exclusively to the region occupied by the O₂-O₂ profile, 465–491 nm. The line edges are further refined by applying piecewise fits (first- or second-degree polynomials) to the relatively narrow windows, 459–465 and 486–491 nm, thus concluding this rather involved procedure for the first category of the low-reflectance scenes.

Step 4, the SCD retrieval, follows the approach described in Marchenko et al. (2015). Here, preliminary SCD values are obtained from two algorithms, the Nelder–Mead minimization method and the least-squares Levenberg–Marquardt fit (Press et al., 1992), taking the latter as a default and fitting the normalized (Step 3) O₂-O₂ profile in the 465–487 nm interval. These evaluations are repeated for each temperature-dependent O₂-O₂ cross section; there are five of them measured by Thalman and Volkamer (2013). Each cross section is fitted to the data, providing an individual root mean square (RMS) value of the fitting residuals. These five RMS values are approximated by a parabolic function. The minimum of the function is used to construct via linear interpolation a synthetic O₂-O₂ profile that is removed from the normalized reflectances, thus leaving us with the residuals that are presumably dominated by instrumental noise. The noise is reduced by an iterative procedure similar to one described in Marchenko et al. (2015). The final SCD evaluation is performed over a slightly broadened wavelength range, 463–488 nm.

As implemented, the algorithm relies on optimal SCD retrievals of the O₃, NO₂, and H₂O trace gases, as well as preliminary cloud-fraction estimates. The latter is used exclusively over deep-water areas during the wavelength calibration and the Raman scattering removal. If needed, such cloud fractions can be substituted for appropriately adjusted reflectances, thus vying for self-sufficiency. The use of independent O₃, NO₂, and H₂O SCDs is an essential part of the algorithm that, especially for the scenes with heavy O₃ and NO₂ loads, leads to more accurate O₂-O₂ SCDs. The use of the trace-gas SCDs does not create any paradox when the NO₂ values are used in order to retrieve cloud properties that should be incorporated into the NO₂ estimates. Note that in the implemented algorithm we use the NO₂ SCD estimates that can be obtained without any relevance to cloud properties. These cloud properties are used later, during the conversion of the NO₂ slant columns to the NO₂ vertical columns. Opting for a complete self-reliance of the cloud algorithm, one may substitute the required O₃, NO₂, and H₂O SCDs for SCD estimates provided by the appropriate trace-gas climatologies.

Figure 6 ECFs retrieved with our algorithm versus those retrieved from OMCLDO2. Results are provided as 2-D densities in ECF bins of 0.01. The color scale represents the number of OMI pixels falling within a given bin. Data for 13 November 2006, 30^∘ S–30^∘ N. (a) Land; (b) ocean.

Figure 7 Similar to Fig. 6 but for OCP.

Figure 8 Comparison of scene pressures from our algorithm (black curve) and OMCLDO2 v2 (red) with surface pressures (blue) along cross-track position 20 of OMI orbit 4415 (14 May 2005). The green curve shows cloud fractions for this cross-track position.

2.3 Cloud algorithm

The O₂-O₂ cloud algorithm described here is based on the O₂-O₂ absorption band at 477 nm. This algorithm is broadly similar to the operational O₂-O₂ cloud algorithm developed at the Royal Netherlands Meteorological Institute (KNMI) known as OMCLDO2 (Acarreta et al., 2004; Sneep et al., 2008; Veefkind et al., 2016). However, our approach differs in a number of aspects.

First, we use normalized radiance at 466 nm to compute f with Eq. (1) in a separate step. This wavelength was selected because it is not significantly affected by gaseous absorption and rotational Raman scattering and it is still sufficiently close to the O₂-O₂ absorption band center at 477 nm. f is calculated using linear interpolation of lookup tables (LUTs) of I_g and I_c. The tables were generated for 23 different surface and cloud pressures, 20 surface reflectivities, 30 SZAs, 20 VZAs, and 20 relative azimuth angles. Nodes and their locations were selected on the basis of the analysis of interpolation errors. A threshold for acceptable interpolation error was set at 0.2 %. It should be noted that aerosols are implicitly accounted for in the determination of f, as they are treated (like clouds) as particulate scatters.

Figure 9 Difference between the scene pressure and the surface pressure, P_sc−P_s, for our algorithm (a) and OMCLDO2 v2 (b). Data are for 13 November 2006.

Figure 10 Two-dimensional histograms comparing effective cloud fraction (ECF) retrieved with GLER (y axes) and climatological LER (x axes) for land (a) and ocean (b). The color scale shows numbers of data points. OMI data are for 13 November 2006.

Figure 11 Differences between ECFs retrieved with GLER and those retrieved with climatological LER, f(GLER)−f(ClimLER). OMI data are for 13 November 2006. No snow-/ice-covered areas are included in the comparison.

Figure 12 Cross-track dependence of ECF zonal means retrieved with GLER (black) and climatological LER (red) for different latitude bins. Data are for 13 November 2006; effective cloud fractions are between 0.05 and 0.25.

Our algorithm retrieves cloud OCP from the OMI-derived oxygen dimer SCD at 477 nm. The OCP, here also denoted as P_c, is estimated using the MLER method to compute the appropriate air mass factors (AMFs) (Yang et al., 2015). To solve for OCP, we invert the following equation:

$$\begin{array}{ll}{\displaystyle }\mathrm{SCD}& {\displaystyle }={\mathrm{AMF}}_{\mathrm{g}}(P_{\mathrm{s}},R_{\mathrm{g}})\mathrm{VCD}\left(P_{\mathrm{s}}\right)(\mathrm{1}-f_{\mathrm{r}})\\ \text{(4)}& {\displaystyle }& {\displaystyle }+{\mathrm{AMF}}_{\mathrm{c}}(P_{\mathrm{c}},R_{\mathrm{c}})\mathrm{VCD}\left(P_{\mathrm{c}}\right)f_{\mathrm{r}},\end{array}$$

where VCD is the vertical column density of O₂-O₂ ($\mathrm{VCD}=\mathrm{SCD}/\mathrm{AMF}$), AMF_g and AMF_c are the precomputed (at 477 nm) clear-sky (subscript g) and overcast (cloudy, subscript c) subpixel AMFs, P_s is the surface pressure, and f_r is the cloud radiance fraction (CRF) given by $f_{\mathrm{r}}=f{I}_{\mathrm{c}}/I_{\mathrm{m}}$. Equation (4) is similar to that frequently used for retrieval of trace-gas VCDs with the MLER model provided P_c and R_c are known (see, e.g., Veefkind et al., 2006). Here we use this equation for retrieval of P_c assuming that the O₂-O₂ VCD is known. The CRF is calculated at 466 nm. CRF defines a fraction of TOA radiance reflected by the cloud. It should be noted that CRF is wavelength dependent (see discussion in Sect. 3.4.1). The CRF retrievals at different wavelengths are included in our output. Lookup tables of the TOA radiances and AMFs were generated using VLIDORT. Temperature profiles needed for estimation of VCD and AMF are taken from the NASA GMAO GEOS-5 global data assimilation system (Rienecker et al., 2011).

In addition to OCP, we retrieve the so-called scene pressure, P_sc. The scene pressure is derived from Eq. (4) assuming that f_r=1 and R_c is equal to the scene LER, R_sc:

$$\begin{array}{}\text{(5)}& \mathrm{SCD}={\mathrm{AMF}}_{\mathrm{c}}(P_{\mathrm{sc}},R_{\mathrm{sc}})\mathrm{VCD}\left(P_{\mathrm{sc}}\right).\end{array}$$

R_sc is determined from the measured TOA radiance using Eq. (2) for a known surface pressure. In the absence of clouds and aerosols, the P_sc should be equal to P_s. P_sc is therefore an important diagnostic tool for evaluation of the performance of cloud pressure algorithms.

Figure 13 Two-dimensional histogram similar to Fig. 10 but comparing OCPs retrieved with GLER with those retrieved with climatological LER.

Figure 14 Differences between OCPs retrieved with GLER and climatological LER for 13 November 2006. Data are shown for ECF>0.05 only.

Figure 15 Cross-track dependence of OCPs retrieved with GLER (black) and those retrieved with climatological LER (red) for 13 November 2006.

3 Results and discussion

3.1 Evaluation of the cloud algorithm

To evaluate our cloud algorithm we have compared the retrieved values of f and P_c with those from the operational OMCLDO2 version 2 (Veefkind et al., 2016). For this comparison, the cloud products are retrieved for the climatological surface LER (Kleipool et al., 2008) identical to that used in OMCLDO2 v2. Figure 6 shows scatter plots of effective cloud fractions (ECFs) calculated with our algorithm versus those calculated with OMCLDO2 for land and ocean for a selected date, 13 November 2006. Different dates show very similar trends. The scatter of data around the 1:1 line is somewhat higher for low values of fs. The mean differences between the two data sets do not exceed 0.01 for all values of f. The standard deviation of the f differences is within 0.01 for ocean and 0.03 for land. Differences in values of f are probably due to contrasting approaches used in the two algorithms. We retrieve the ECF at 466 nm independently from the OCP retrieval, whereas OMCLDO2 retrieves the ECF and OCP simultaneously at 477 nm.

Figure 7 shows scatter plots of OCPs calculated with our algorithm versus those calculated with OMCLDO2 v2 for the same day. There is a bias between OMCLDO2 and our algorithm: our OCP retrievals are higher than those from OMCLDO2 by about 50 hPa on average. The standard deviation of the OCP differences ranges from about 100 hPa for OCP<400 to 150 hPa for lower OCPs. Higher OCP retrievals from our algorithm compared to OMCLDO2 can be related to slightly higher O₂–O₂ SCD estimates and also to differences in ECF, which affect the OCP retrievals.

Figure 8 shows a comparison of P_s and P_sc from our algorithm and OMCLDO2 v2 along cross-track position 20 of the OMI orbit 4415. A similar comparison of the OMCLDO2 v2 P_sc and with P_s for this cross-track position was carried out in Veefkind et al. (2016). We added our scene pressure and cloud fraction for this cross-track position to this comparison. First, P_sc retrievals from both algorithms agree very well with the surface pressures for the high-reflectivity scenes in Greenland (OMI scan lines along the orbit with numbers iTimes=1300–1400). For mostly cloud-free conditions over the ocean (iTimes∼400 and iTimes ∼270), the scene pressures retrieved from both algorithms are higher than the surface pressures. These differences for those oceanic regions are slightly lower for our algorithm than for OMCLDO2. Over desert and semidesert areas, P_sc retrieved with our algorithm over mostly cloud-free conditions (iTimes=720–900 and iTimes∼570) is close to the surface pressure while OMCLDO2 significantly overestimates the scene pressure. We attribute the better performance of our algorithm over deserts and semideserts to the special adjustment of the spectral fitting procedure for those areas (see Sect. 2.2 and Fig. 5). Figure 9 shows maps of differences between the scene pressure and the surface pressure retrieved from our algorithm and OMCLDO2 v2. The differences are shown for mostly clear scenes with f<0.25. A comparison of the maps shows that the performance of our algorithm over land is better than the performance of OMCLDO2 v2. This is clearly seen for the Sahara, the Arabian Peninsula, Australia, etc. Over the ocean, our algorithm performs slightly better only for certain areas in the Southern Ocean, e.g., for three locations at ∼55^∘ S, ∼60^∘ W, ∼20^∘ W, and ∼100^∘ E. Over other areas of the ocean, the performance of the two algorithms is similar.

3.2 Comparison of ECFs derived with GLER and climatological LER

Figure 10 shows scatter plots of f calculated with GLER versus those calculated with climatological LER (ClimLER) (Kleipool et al., 2008) in the form of 2-D histograms. The color scale on the 2-D histograms represents a number of data points. As expected, the scatter and systematic deviation of data around the 1:1 line diminish with increasing f. Over land, f calculated with GLER is mostly higher than f calculated with the climatological LER (Fig. 10a) particularly for low f. This is explained by differences between GLERs and climatological LERs (Vasilkov et al., 2017), which have the most pronounced impact on f for low cloudiness. The GLER values are mostly lower than the climatological LERs because the former are derived from atmospherically corrected MODIS radiances while the latter are affected by residual aerosols. Moreover, climatological LERs can be contaminated by clouds owing to substantially larger sizes of OMI pixels compared with those of MODIS. As it follows from Eq. (1), lower surface LER leads to lower clear-sky radiance, thus increasing f. For the most important NO₂ retrieval range of f<0.25, f retrieved with GLER is higher than that retrieved with climatological LER by approximately 0.02 on average. The standard deviation of the GLER–ClimLER f differences varies between 0.03 and 0.05 depending on the f value.

Over ocean, the GLERs are higher than the climatological LERs in the areas affected by sun glint and areas observed at large viewing zenith angles (Vasilkov et al., 2017). Therefore, f retrieved with GLER in these areas is lower than that retrieved with the climatological LER (Fig. 10b). For other oceanic regions, the GLERs are slightly lower than the climatological LERs, resulting in slightly higher values of f. For the range of f<0.25, which is frequently used in tropospheric trace-gas retrievals, the mean difference between f retrieved with GLER and that retrieved with climatological LER is approximately 0.02 when averaged globally. The standard deviation of the f differences varies within 0.02–0.03 depending on the f value. Even though the f differences are small on average, they can be as large as 0.05–0.07 for individual pixels – this is quite significant for the low f range. It should be noted that the fraction of negative f retrievals is lower when using GLER compared with the climatological LER. This is clear evidence of the improvement of f retrievals with GLER.

Figure 16 Differences between ECFs retrieved at 405 and 435 nm and ECF retrieved at 466 nm as a function of ECF at 466 nm for 13 November 2006. (a) Land; (b) ocean.

Figure 17 Similar to Fig. 16 but for CRF.

Figure 18 (a) Monthly mean gridded ($\mathrm{0.5}{}^{\circ }\mathrm{latitude}\times \mathrm{0.5}{}^{\circ }$ longitude) OMI NO₂ tropospheric VCDs for January 2006 retrieved using the GLER and O₂-O₂ cloud products. (b) Change in tropospheric VCDs due to the change in surface reflectivity for NO₂ retrievals alone. (c) Change in tropospheric VCDs due to the change in surface reflectivity for both NO₂ and cloud retrievals.

Figure 19 Similar to Fig. 18 but for July 2006.

Figure 11 is a geographic map of differences between f calculated with GLER and climatological LER for 13 November 2006. Over the ocean, the most prominent features are sun glint areas where the negative f differences (GLER–ClimLER) are at maximum. This is because the climatological LERs are derived from minimum values of LERs from a long time series (up to 5 years) of observations over a given area; that is why the impact on observations affected by sun glint is somewhat mitigated. Over land, the f differences are mostly positive due to aerosol and possible cloud contamination of the climatological LERs (note that the relatively large OMI footprints lead to predominance of cloudy scenes).

One indicator of the cloud algorithm performance is the cross-track dependence of the retrieved f and OCP averaged over a given latitude bin. Ideally, this dependence should be fairly flat if the retrieved parameters are averaged globally. Figure 12 compares the cross-track dependencies of f retrieved with GLER and ClimLER for 13 November 2006. The comparison is carried out for 30^∘ latitude bins for land and ocean separately and specifically for the low f range $\mathrm{0.05}<f<\mathrm{0.25}$. Over land, the cross-track dependence of f is reasonably flat. f retrieved with GLER is mostly lower than that retrieved with the climatological LERs because GLER is generally lower than the climatological LER. Over the ocean, the most noticeable feature of the cross-track dependence is related to the sun glint area in the 30–0^∘ S latitude bin. Even though the use of GLER smooths out the irregularity in the cross-track dependence of f, it slightly overcorrects this irregularity owing to overestimation of GLER in sun glint areas. We have not yet determined the exact cause of this GLER overestimation. However, our preliminary radiative transfer simulations show that the presence of nonabsorbent aerosol with low optical thickness (τ<0.2) can reduce the sun glint GLER. It should be noted that the cross-track dependence of f for cloudy scenes with f>0.25 is much flatter than that for scenes with f<0.25.

3.3 Comparison of OCPs produced with GLER and climatological LER

Figure 13 shows 2-D histogram plots of OCPs calculated with GLER and climatological LER. Over land, the GLER-retrieved OCPs are slightly higher than those retrieved with climatological LER (Fig. 13a). Over ocean, the data are closer to the 1:1 line (Fig. 13b) but the overall scatter is somewhat higher than that over land. Data for which the GLER-retrieved OCPs are lower than those retrieved with climatological LER are mostly from sun glint areas. There is a cluster of data for which both GLER-retrieved OCPs and those retrieved with climatological LER are higher than the surface pressure (1013 hPa). These data are retrieved over virtually clear-sky conditions, and there are two possible reasons for this effect. First, OCP retrievals higher than the surface pressure are evidence of enhanced O₂-O₂ absorption. This can be caused by scattering, low-altitude aerosols in which enhanced photon path length prevails over other aerosol effects. A second possible reason is due to a remaining deficiency in our spectral fitting procedure, which somehow overestimates SCDs over the clear-sky ocean. A very small fraction of the OCP retrievals with values around 100 hPa are likely due to an artifact caused by the LUT extrapolation over the minimum pressure node of 100 hPa.

Geographic distribution of the differences between OCPs calculated with GLER and climatological LER is shown in Fig. 14. Over ocean, most areas with negative differences (OCP(GLER) – OCP(ClimLER)) correspond to sun glint.

Figure 15 shows cross-track dependencies of OCPs retrieved with GLER and climatological LER for 13 November 2006. Again, the range $\mathrm{0.05}<f<\mathrm{0.25}$ is used for the analysis. Over the ocean areas affected by sun glint, OCPs retrieved with GLER are significantly lower than those retrieved with climatological LER. The underestimation of OCP is mainly related to underestimation of O₂-O₂ SCDs in those areas. This issue requires further investigation. It should be noted that the cross-track dependence of OCP for cloudier scenes with f>0.25 is much flatter than that for lower cloudy scenes with f<0.25. The cross-track dependence of OCPs significantly deviates from the uniform dependence for pixels near the edges of the OMI swath. The OCPs at the edges are substantially higher than those in the near-nadir parts of the swath. This behavior of the OCP cross-track dependence is observed for both land and ocean and is not affected by the use of different surface reflectivities. This increase in OCP retrievals at the swath edges is not understood. It should be noted that similar behavior of the OCP cross-track dependence is also seen in the operational OMCLDO2 v2 retrievals for low f.

3.4 Application to the OMI NO₂ algorithm

3.4.1 Spectral dependence of ECF and CRF

The OMI NO₂ algorithm uses the DOAS approach to fit OMI-measured spectra in the wide spectral window of 405–465 nm. Hence, the question arises as to how the value of f determined at 466 nm can be representative for the entire spectral window. A simulation experiment carried out by Gupta et al. (2016) showed that the wavelength dependence of ECF is weak. In this experiment, observed TOA radiances were simulated as a sum of the clear-sky and cloudy radiances weighted with a geometrical cloud fraction. It was shown that f varied only a few percent over a wide spectral range from UV to near IR. We have verified this result using OMI-observed spectra and the calculation of f at two additional wavelengths: 405 and 435 nm. A lookup table of TOA radiances at these wavelengths was generated using VLIDORT. To calculate the clear subpixel TOA radiance we use the climatological surface reflectivity from Kleipool et al. (2008) with linear interpolation of the spectral dependence of the surface reflectivity.

Figure 20 Difference in tropospheric NO₂ VCD caused by replacing climatological LER with GLER for July (a) and January (b) 2006. Surface reflectivity affects NO₂ retrievals directly as an input to the AMF calculation (blue symbols) and indirectly by changing cloud parameters used in the AMF calculation (red symbols). Vertical bars represent the standard deviation for each class of tropospheric NO₂ VCD of the size 1×10¹⁵ molec cm⁻². The standard deviations contain both effects of surface reflectivity on NO₂ retrievals.

Figure 16 shows differences between the baseline f at 466 nm and f calculated at additional wavelengths along with the standard deviation. The differences are shown as a function of f for land and ocean separately. The difference of f calculated at 435 nm is obviously lower than that for f calculated at 405 nm. The f differences decrease with increasing f due to the gradually decreasing contribution of the clear subpixel to the TOA radiance, which is mostly responsible for the f spectral dependence. The f differences are at maximum (less than 0.02) in the middle range of $\mathrm{0.4}<f<\mathrm{0.6}$. The f differences decrease with decreasing f. However, the relative f differences increase with decreasing f because of decreasing absolute values of f. Over the ocean, the f differences are slightly higher than those for land except for low values of f. For f<0.1 the f differences over ocean are noticeably lower than those over land. Overall, the spectral differences of f within the NO₂ retrieval window are small and do not exceed ∼0.01 over land and ∼0.015 over ocean for the most important range of f<0.25.

CRF is used for calculation of trace-gas AMFs in cloudy conditions. The CRF dependence on wavelength is much more pronounced than the spectral dependence of f mostly due to the spectral dependence of the measured TOA radiances. For thick cloud, CRF is simply a fraction of TOA radiance reflected by the cloud. The CRF varies with wavelength because the radiance coming from the cloud-free part of the scene is wavelength dependent. A physical interpretation of the CRF for thin clouds is not obvious. Figure 17 shows differences between CRFs calculated at 466 nm and CRFs calculated at 405 and 435 nm. Over the ocean, the f differences between 466 and 405 nm are somewhat higher than those over land. For f<0.25, the CRF differences do not exceed ∼0.07–0.08 within the NO₂ retrieval window. Such differences are quite acceptable for purposes of tropospheric NO₂ retrievals.

3.4.2 Cloud effects on NO₂ retrievals

We used the operational OMI NO₂ algorithm, OMNO2 version 3 (Krotkov et al., 2017), to assess how the change in surface reflectivity affects the retrievals of stratospheric and tropospheric NO₂ VCDs. This algorithm comprises four main steps: (1) retrieval of NO₂ SCDs by spectral fitting of laboratory-measured spectra to the OMI-measured absorption spectrum in the range of 402–465 nm (Marchenko et al., 2015); (2) calculation of AMFs using various input parameters such as viewing geometry, surface reflectivity, cloud pressure, cloud radiance fraction, and a priori NO₂ profile shapes; (3) removal of cross-track striping; and (4) conversion of SCDs to VCDs using AMFs and separation of stratospheric and tropospheric components (Bucsela et al., 2013). Since the retrieval of cloud parameters (cloud pressure and ECF) is also affected by surface reflectivity, changes in surface reflectivity affect NO₂ retrievals both directly, as inputs to the AMF calculation, and indirectly, through the cloud parameters.

We conducted separate NO₂ retrievals using climatological LER (Kleipool et al., 2008) and GLER (Vasilkov et al., 2017) and cloud parameters retrieved using the respective surface reflectivity products. Analysis of stratospheric NO₂ VCDs revealed that changing surface reflectivity only in NO₂ retrievals (i.e., no change in cloud parameters) had a minor impact (<1 %) on stratospheric NO₂ estimates. The effect was somewhat larger for the changes in cloud parameters, with the difference in estimated stratospheric NO₂ VCDs reaching up to 5 %. This is expected because the troposphere–stratosphere separation scheme in OMNO2 uses NO₂ observations from unpolluted and cloudy areas to construct the stratospheric NO₂ field.

Figures 18a and 19a show monthly mean GLER-based tropospheric NO₂ VCDs for July and January 2006, respectively. Tropospheric NO₂ exhibits strong spatial variability, with pronounced enhancements over industrial and other source regions. In addition, we observe higher tropospheric NO₂ in January compared to July, a reflection of the relatively longer NO₂ lifetime and shallower boundary layer in winter. In contrast to stratospheric NO₂ VCDs, retrievals of tropospheric NO₂ VCDs are very sensitive to the changes in surface reflectivity and cloud parameters. The bottom two panels of Fig. 18 and 19 show how retrievals of tropospheric NO₂ VCDs are affected by replacing surface reflectivity from climatological LER with GLER. Using GLER-based surface reflectivity reduces tropospheric AMFs and enhances tropospheric NO₂ VCDs. The impact of changes in surface reflectivity varies with the vertical distribution of NO₂, with the largest effects in polluted areas.

Figure 20 quantifies the percent change in tropospheric NO₂ VCDs as a function of NO₂ levels, suggesting the GLER effect on NO₂ retrievals can reach as much as 20–30 %. Additional effects of LER changes on NO₂ retrievals come through changes in cloud parameters that may cause 10–15 % of additional changes in tropospheric NO₂ VCDs.

4 Conclusions

We have developed a new cloud algorithm based on the O₂-O₂ absorption band at 477 nm. The main features of the algorithm are (1) a new spectral fitting method of TOA radiances to derive O₂-O₂ SCDs and (2) the use of surface GLERs that replaces climatological LER data sets. This new cloud algorithm is intended for use within the standard OMI NO₂ algorithm and planned as a backup algorithm for the upcoming TEMPO geostationary mission.

Validation of our cloud algorithm was carried out by comparisons of the retrieved values of f and OCP with values from the latest version of the OMI operational algorithm, OMCLDO2 v2, also based on the O₂-O₂ absorption band at 477 nm. f and OCP were retrieved for the climatological surface LER identical to OMCLDO2. Comparisons showed a good agreement between our f and that from OMCLDO2. Our OCPs are overall higher than those from OMCLDO2 by about ∼50 hPa on average. Diagnostic scene pressures from our algorithm are slightly closer to the surface pressure than those from OMCLDO2.

We examined f and OCP changes caused by replacing the climatological surface LERs with GLERs. For the scenes with f<0.25, the range traditionally used in the trace-gas retrievals, values of f retrieved with GLER are higher than those retrieved with climatological LER by ∼0.02 on average. Even though the f differences are small on average, they can be as large as 0.05–0.07 for individual pixels; this is quite significant for the low values of f. Over land, the GLER-retrieved OCPs are slightly higher than those retrieved with climatological LER. Over ocean, the data are closer to the 1:1 line than those over land. The geographical regions where the GLER-retrieved OCPs are lower than those retrieved with climatological LER are mostly related to sun glint areas.

We applied the new cloud algorithm to OMI NO₂ retrievals and analyzed NO₂ column changes related to the use of the new cloud algorithm. The GLER effect on NO₂ AMFs can increase tropospheric NO₂ retrievals by 20–30 % over polluted regions. An effect on NO₂ retrievals that comes through changes in cloud fraction and pressure can make 10–15 % of additional changes in tropospheric NO₂ VCDs.