Articles | Volume 15, issue 20
Research article
19 Oct 2022
Research article |  | 19 Oct 2022

Algorithm theoretical basis for ozone and sulfur dioxide retrievals from DSCOVR EPIC

Xinzhou Huang and Kai Yang

On board the Deep Space Climate Observatory (DSCOVR), the first Earth-observing satellite at the L1 point (the first Lagrangian point in the Earth–Sun system), the Earth Polychromatic Imaging Camera (EPIC) continuously observes the entire sunlit face of the Earth. EPIC measures the solar backscattered and reflected radiances in 10 discrete spectral channels, four of which are in the ultraviolet (UV) range. These UV bands are selected primarily for total ozone (O3) and aerosol retrievals based on heritage algorithms developed for the series of Total Ozone Mapping Spectrometers (TOMS). These UV measurements also provide sensitive detection of sulfur dioxide (SO2) and volcanic ash, both of which may be episodically injected into the atmosphere during explosive volcanic eruptions. This paper presents the theoretical basis and mathematical procedures for the direct vertical column fitting (DVCF) algorithm used for retrieving total vertical columns of O3 and SO2 from DSCOVR EPIC. This paper describes algorithm advances, including an improved O3 profile representation that enables profile adjustments from multiple spectral measurements and the spatial optimal estimation (SOE) scheme that reduces O3 artifacts resulting from EPIC's band-to-band misregistrations. Furthermore, this paper discusses detailed error analyses and presents intercomparisons with correlative data to validate O3 and SO2 retrievals from EPIC.

1 Introduction

The Deep Space Climate Observatory (DSCOVR) was launched on 11 February 2015 and after a 116 d journey successfully maneuvered into its Lissajous orbit around the first Earth–Sun system Lagrangian (L1) point, which is about 1.5×106 km from the Earth and located between the Sun and the Earth on the ecliptic plane. At the L1 point, where the net gravitational forces equal the centrifugal force, DSCOVR orbits the Sun at the same rate as the Earth, staying closely in line along the Sun and the Earth, thus allowing the Earth-pointing EPIC to continuously monitor the entire sunlit planet.

The Earth Polychromatic Imaging Camera (EPIC) measures the solar backscattered and reflected radiances from the Earth using a two-dimensional (2048×2048) charged-coupled device (CCD) recording a set of 10 spectral images successively using different narrowband filters. While EPIC may continuously observe the Earth from the vicinity of the L1 point, only a number of spectral image sets are taken in a day, limited by accessible contact windows of the two ground stations located in Wallops island (Virginia) and Fairbanks (Alaska). Currently, between 13 and 22 spectral image sets, recorded at a sampling rate of one set in every 110 min during boreal winter and every 65 min during boreal summer, are transmitted back to the ground stations in a day.

EPIC takes about 6.5 min to complete an image set. The first in the set is the blue band (centered at 443 nm), which takes ∼2 min to complete the imaging at native resolution (2048×2048 pixels). The images of the nine remaining bands are sequentially recorded at a reduced resolution (1024×1024 pixels, achieved through an onboard average of 2×2 pixels), separated by a time cadence of ∼30 s between adjacent bands. Due to the Earth's rotation and spacecraft jitter, each spectral image records a slightly different (i.e., rotated) sunlit hemisphere. As a result, the images of two different channels appear to be displaced from each other, usually by a distance of about one to a few native pixels, depending on their observation time difference.

Each native pixel has a ∼1 arcsec or 2.778×10-4 angular instantaneous field of view (IFOV), yielding a geometric ground footprint size of 8×8 km2 at the image center of the sunlit disk. The effective footprint size is about 10×10 km2, which is larger than the geometric one due to the effect of the optical point-spread function of the EPIC imaging system. For a reduced-resolution image (1024×1024 pixels), the effective central ground IFOV size is about 18×18 km2, which is significantly smaller than the nadir footprints of some past and present satellite instruments that provided global ozone mapping from the low Earth orbit (LEO), such as the Total Ozone Mapping Spectrometer (TOMS; nadir pixel size 50×50 km2) on a series of satellites, the Scanning Imaging Absorption Spectrometer for Atmospheric Cartography (SCIAMACHY; 60×30 km2; Bovensmann et al.1999) on ESA's ENVIronmental SATellite (ENVISAT), the Ozone Mapping and Profiler Suite Nadir Mapper (OMPS-NM; 50×50 km2; Flynn et al.2014) on the Suomi National Polar Partnership (SNPP), and the Global Ozone Monitoring Experiment-2 (GOME-2; Callies et al.2000; Munro et al.2016) on Metop-A (40×40 km2), Metop-B (80×40 km2), and Metop-C (80×40 km2). Though it is slightly larger than the nadir footprint of the Ozone Monitoring Instrument (OMI; 13×24 km2; Levelt et al.2006) on Aura and the OMPS-NM (17×13 km2; Flynn et al.2016) on NOAA-20, as well as much bigger than that of the TROPOspheric Monitoring Instrument (TROPOMI; 5.5×3.5 km2; Veefkind et al.2012) on the ESA Sentinel-5 Precursor (S5P), EPIC's spatial resolution is sufficiently high to map small-scale O3 natural variations and observe many volcanic emissions, from degassing to eruption.

EPIC, combining moderate spatial resolution with high temporal cadences from the unique vantage point of L1, provides unprecedented Earth observations from sunrise to sunset simultaneously (see Fig. 1). This synoptic (i.e., concurrent, globally unified, and spatially resolved) perspective is quite distinctive from satellite observations from an LEO or geostationary Earth orbit (GEO): LEO observations are often made within a narrow range of local time with a small number of samplings at a location per day, while GEO observations have limited spatial coverage, constrained to roughly 60 away from its position. The EPIC observations can have simultaneous co-located observations with measurements from any contemporaneous LEO and GEO platforms, allowing direct comparisons and synergistic use of data acquired from different perspectives. This overlapping feature has been exploited to calibrate some EPIC channels by matching its measured albedo values to those of OMPS-NM on SNPP (Herman et al.2018).

Figure 1(a) Example of the EPIC field of view (FOV): EPIC Earth image at 11:40:31 UTC on 4 September 2015. Image source: NASA EPIC Team, accessed via (last access: 1 October 2022). (b) Viewing and illumination angles are taken from FOV on the left. The subsolar point is marked on the map with a yellow dot. The area shaded with midnight blue is in the dark, i.e., without direct sunlight, while the unshaded area is the sunlit hemisphere, with sunrise on the left (west of subsolar point) and sunset on the right (east of subsolar point). Contours of solar zenith angles (SZAs, blue dashed lines) and viewing zenith angles (VZAs, red dashed lines), going from 10 to 80 with a step of 10, are shown in the sunlit area. Note that the SZA (θs) and VZA (θv) of an EPIC IFOV have similar values, and both angles increase as the IFOV moves from the center towards the edge of the sunlit disk.

The 10 narrow bands of EPIC, spanning ultraviolet (UV), visible, and near-infrared wavelengths, are selected to yield diverse information about the Earth, from atmospheric compositions to surface reflectivity and vegetation. Four of the 10 bands measure UV spectral radiances, which are primarily used for total ozone (O3) retrievals. These UV bands also provide sensitive detection of sulfur dioxide (SO2) and volcanic ash, both of which may be episodically injected into the atmosphere during explosive volcanic eruptions.

This paper describes algorithm physics, model assumptions, mathematical procedures, and error analyses for the direct vertical fitting (DVCF) algorithm. We show examples to illustrate the high accuracy of O3 and SO2 retrievals achieved by applying the DVCF algorithm to spectral UV radiance measurements of DSCOVR EPIC. Lastly, we validate the DSCOVR EPIC O3 and SO2 through intercomparisons with correlative data.

2 Algorithm physics

Algorithm physics is a term first used by Chance (2006) to denote the physical processes contributing to the spaceborne measurement of radiance spectra. A measured radiance Lm (in units of Wsr-1m-2nm-1) from space consists of sunlight photons within a narrow spectral range (typically <2 nm), specified by the instrument spectral response function S (ISRF, e.g., EPIC UV filter transmissions shown in Fig. 2), and is modeled as

(1) L M = S ( λ ) I TOA ( λ ) F ( λ ) d λ S ( λ ) d λ ,

where F(λ) (in units of Wm-2nm-1) is the monochromatic spectral solar irradiance, and ITOA(λ) is the Sun-normalized monochromatic top-of-the-atmosphere (TOA) radiance (in units of sr−1) for a wavelength λ (in units of nm). The Sun-normalized measured radiance IM for a spectral band is defined as IM=LM/FM, where FM=S(λ)F(λ)dλ/S(λ)dλ, and the λ integrations in these equations are performed over the valid range of the ISRF S for the spectral band. Hereafter we drop “Sun-normalized” when referring to IM, which is simply called measured radiance. Quantities for a spectral band are flux-weighted bandpass averages to account for the differential contributions from individual wavelengths within the bandpass. Without loss of generality, ITOA(λ) and other spectral-dependent quantities are hereafter used to denote flux-weighted bandpass averages, with λ representing the characterized wavelength of the spectral band.

Figure 2Filter transmission functions for the four EPIC UV channels. The widths are ∼1 nm for EPIC bands 1 and 2, similar to those for TOMS and OMPS-NM. Note that the filter transmissions as functions of wavelength are measured in the air (see Fig. 1 in Herman et al.2018). Here we have converted the wavelength in the air to wavelength in a vacuum using the formula of Edlén (1966). The filter values are normalized to 1 at band centers (noted on top of each panel with uncertainty).


To reach a sensor at TOA, sunlight photons are either back-scattered by air molecules or particles or reflected by the underlying Earth surface. As these photons traverse through the atmosphere along many possible optical paths connecting the Sun to the sensor, they may be absorbed by the underlying surface or by some atmospheric constituents, such as trace gases (e.g., O3 and SO2) and light-absorbing particles (e.g., dust and smoke). The photons that complete the journey carry information about atmospheric absorbers along their paths. The accumulation of photons from each contributing path yields the TOA radiance, which may be modeled with radiative transfer (RT) simulation if the properties of surface reflection as well as atmospheric absorption and scattering are known explicitly. The ability to model the TOA radiance accurately is the prerequisite for interpreting the observations and relating the gas absorptions with TOA radiance measurements.

We next describe the characteristics of UV photon sampling of the atmosphere and the construction of surface and atmospheric models to enable proper simulation of the photon sampling of the atmosphere. Dividing the atmosphere into infinitesimal thin layers, the quantity that specifies the photon sampling is the mean path length of photons traversing through a layer. This mean path length normalized by the geometric thickness of the layer is the local or altitude-resolved air mass factor (AMF, mz). The proper simulation of photon sampling requires the modeled mean path length through each layer to closely match that in the actual observing condition.

In theory, a TOA radiance, ITOA, depends on the viewing illumination geometry, the optical properties of the atmospheric constituents (both absorbers and non-absorbers), and their amounts and vertical distributions, as well as on the reflective properties of the underlying surface. For a wavelength λ, ITOA can be expressed as the sum of two contributions,

(2) I TOA = I a + I s ,

where Ia consists of solar photons scattered once or more by molecules and particles in the atmosphere without interacting with the underlying surface, and Is represents solar photons reflected at least once or multiple times by the underlying surface.

2.1 Path radiance

Ia is also known as the atmospheric path radiance, i.e., photons backscattered to the sensor along a path without any intersection with the underlying surface. Conceptually it is the accumulation of TOA photons that are last backscattered toward the sensor along the line of sight from atmospheric layers at different levels of extinction optical depths. Algebraically it is expressed as the path integration of virtual emission J(t) (Dave1964) in the direction specified by the view zenith angle (θv), attenuated (e-t/μ, where μ=cos θv) by atmospheric scattering and absorption, over the extinction optical depth t along the path of line of sight from the top (t=0) to the bottom (t=τ) of the atmosphere:

(3) I a = 0 τ J ( t ) e - t / μ ω ( t ) d t / μ .

The source of virtual emission, J(t), consists of all the photons scattered towards to the sensor, including photons of the direct solar radiation being scattered once only and photons of diffuse radiation (i.e., photons scattered to level t) being scattered once more at t. The strength of the virtual emission of a thin layer at t is proportional to its scattering optical thickness, which is equal to the product of the layer total optical thickness (dt) and the single-scattering albedo ω(t) (defined as the ratio of layer scattering optical thickness over the layer total optical thickness). Here we use Ψ(t)=J(t)e-t/μω(t)/μ to represent the radiance contribution per unit optical thickness to Ia from a layer at t. Equation (3) describes how the solar photons sample the atmosphere from top to bottom and how atmospheric absorption is directly imprinted (via the attenuation e-t/μ) on the path radiance.

A path radiance Ia for a molecular (i.e., an aerosol- and cloud-free) atmosphere with absorption from trace gases can be accurately determined with RT simulations. For example, the path radiances for the low and high zenith angle geometries (see Fig. 3b) are calculated with a vector RT code (e.g., TOMRAD, Dave1964, or VLIDORT, Spurr2006) as a function of wavelength for a molecular atmosphere with the O3 profile X1 in Fig. 3a, and the corresponding radiance contributions to the path radiances at EPIC bands 1 and 2 are shown in Fig. 3c. The radiance contribution function (RCF) for a wavelength in the UV range (300–400 nm) is determined by Rayleigh scattering and absorption by trace gases (primarily O3). O3 is ubiquitous in the atmosphere, with the bulk of it located in the stratosphere (e.g., Figs. 3a or 11), and its absorption cross-sections σ(O3) increase rapidly with shorter wavelengths in the UV range (see Fig. 13). Rayleigh scattering, whose cross-sections are proportional to 1λ4, also increases with shorter wavelength. The strong O3 absorption and large Rayleigh cross-sections at short wavelengths greatly reduce the number of solar photons reaching the lower atmosphere. Conversely, at longer wavelengths, weaker O3 absorption and smaller Rayleigh cross-sections allow more solar photons to reach the lower atmosphere where higher air density increases the intensity of backscattering. Similar to the effect of reducing wavelength, lengthening the slant path (by increasing solar, viewing, or both zenith angles) would enhance ozone absorption and Rayleigh scattering along the slant path, raising the altitude profile of RCF. These spectral and angular characteristics of RCF are illustrated in Fig. 3c, which shows the normalized RCFs (ψ=Ψ/Ia) of EPIC bands 1 and 2 for two different observation geometries and a midlatitude O3 profile labeled as X1 in Fig. 3a. The results in Fig. 3c show that at longer wavelengths and lower zenith angles, path radiance contains more photons that are backscattered from the lower atmosphere. The RCF peak reaches ∼4 km altitude for band 2 at 5 zenith angle, while at shorter wavelength and higher zenith angle, the RCF peak moves to the higher altitude, and it rises to ∼10 km for band 1 at 70 zenith angle. The shifting shapes of RCF shown in Fig. 3c illustrate the changes in the photon sampling of the atmosphere with different wavelengths and zenith angles. The rising RCF peak position signifies diminishing sensitivity to absorptions below the peak while favoring those above it.

Figure 3Sample results from RT simulations for a molecular atmosphere with O3 profiles X1 and X2 in panel (a). Both X1 and X2 are midlatitude zone (30 latitude 60) climatological O3 profiles with the same total vertical column of 275 Dobson units, where 1 DU =2.69×1016 molec. cm−2. RT simulations are performed for two viewing illumination geometries: (1) low zenith angles of θs=θv=5 as well as a relative azimuthal angle of (RAA), ϕ=45 and (2) high zenith angles of θs=θv=70 and ϕ=45. (b) Path radiances Ia(X1) for the low and high zenith geometries and their fractional changes (ΔIa/Ia) when the O3 profile is changed to X2. (c) Normalized RCFs, ψ, for EPIC bands 1 and 2. Here ψ(t) is converted into ψ(lnP) by the multiplication of factor dt/dlnP. (d) Mean photon path lengths (ma) of EPIC bands 1 and 2 as functions of altitude z for the low and high zenith geometries normalized by the respective geometric air mass factors, mG.


The measurement sensitivity to a thin molecular absorber layer is equal to the product of the absorption cross-sections (σ) and the mean path length (ma) of photons passing through the layer, where ma=-lnIa/τz and τz is the absorption optical depth at the layer center altitude z. Note that the photon path length is equal to the geometric AMF, mG=1/cos(θs)+1/cos(θv), for a plane-parallel atmosphere if there is no scattering. Figure 3d shows the mean optical path lengths of EPIC bands 1 and 2 as a function of altitude for the low and high zenith viewing illumination geometries, showing that ma decreases rapidly as the layer descends nearing the surface due to fewer photons reaching the lower atmosphere, while ma approaches mG as the layer rises towards TOA due to fewer path-altering scatterings resulting from lower air density. In the upper troposphere and lower stratosphere (UTLS), ma of the low zenith geometry usually exceeds mG due to a significant fraction of photons undergoing multiple scattering below and within UTLS, while ma of the high zenith geometry drops continuously from TOA down to the surface in the case when the RCF peak is sufficiently high that fewer multiple scatterings contribute to the path radiance. In general, the mean path length ma is shorter for a wavelength with stronger O3 absorption, which reduces the number of photons reaching the lower atmosphere. The variation of ma with a changing altitude signifies the path radiance dependence on the absorber profile. The path radiance fractional change due to profile change (ΔX=X2-X1) can be expressed as

(4) Δ I a I a = I a ( X 2 ) - I a ( X 1 ) I a ( X 1 ) = - 0 Δ X ( z ) σ ( T z ) m a d z ,

where Tz is the atmospheric temperature, and X1(z) and X2(z) are absorber concentration at altitude z. Figure 3b illustrates the change in path radiance caused by an O3 profile change while keeping its total vertical column the same: lowering the O3 profile (e.g., X1 to X2 in Fig. 3a) tends to increase the path radiance. Path radiance changes more with shorter wavelengths at higher zenith angles, thus becoming more sensitive to the shape of the O3 profile. At low zenith angles, the change may have the opposite sign of the change at large zenith angle for certain wavelengths (e.g., the changes plotted as red solid lines for λ>316 nm in Fig. 3b), but the magnitude of change is much smaller, indicating that the path radiances under these conditions are primarily functions of total columns, since they are less sensitive to the profile shapes. The differential responses of the spectral path radiance to profile changes imply that more than one piece of information about O3 may be contained in the multi-spectral measurements. Retrieval constrained by multi-spectral radiances instead of a single spectral band may achieve a more accurate O3 measurement.

2.2 Surface reflection

The path radiance Ia includes backscattered photons that are independent of the underlying surface, while the surface contribution to TOA radiance, Is (referred to as surface radiance hereafter), consists of photons reflected once or more from the surface. For a molecular atmosphere bounded by a surface with well-characterized optical reflection properties, the surface radiance Is can be accurately predicted with RT modeling. For a Lambertian surface, which reflects radiation isotropically independent of the incident direction, the surface radiance Is can be expressed as (Dave1964)

(5) I s = T r s T 1 - r s S b ,

where rs is the reflectance or albedo of the Lambertian surface, T is the total (direct and diffuse) transmittance from the Sun to the surface along the direction of incoming solar irradiation and T from the surface to the TOA along the viewing direction, and Sb is the atmospheric spherical albedo, which is the fraction of the reflected radiation backscattered from the overlaying atmosphere to the surface. The surface contribution from the Lambertian surface, Is, may be described as the once-reflected radiance (TrsT), enhanced by the series of interactions: backscattering from the overlaying atmosphere and reflection from the underlying surface, which are accumulated to produce the amplification factor 1/(1-rsSb).

The reflection property of a surface is represented by a bidirectional reflectance distribution function (BRDF), which specifies the angular distribution of reflected radiance as a fraction of directional incident spectral irradiance. Field measurements (Brennan and Bandeen1970) demonstrate that reflection from natural surfaces (such as cloud, water, and land surfaces) is anisotropic in the UV, exhibiting different apparent reflectances when viewed from different directions. For instance, a water surface looks bright when viewed from the direction near the specular reflection but is much darker outside the glitter (see, e.g., Fig. 4a). Here the apparent reflectance is the Lambertian equivalent reflectivity (LER), i.e., the isotropic reflectance rs that reproduces the radiance Is from a surface with an anisotropic BRDF at a viewing illumination geometry. This LER is also referred to as the geometry-dependent surface LER (GLER) to indicate its dependence on the viewing illumination geometry.

Figure 4Apparent reflectances of an ocean surface described by a Cox–Munk BRDF (Cox and Munk1954a, b) for a wind speed of 6 m s−1, viewed along the plane of incidence with the Sun at a zenith angle of θs=15. (a) GLER at four EPIC UV bands vs. viewing zenith angle θv. Here positive θv denotes ϕ=0 and negative θv for ϕ=180. (b) GLER at several viewing zenith angles vs. wavelength λ.


Reflection of UV sunlight from natural surfaces has long been measured by instruments on board satellites in Sun-synchronous polar orbits (e.g., Eck et al.1987). Since BRDFs for most natural surfaces (except for water surfaces) have not been adequately characterized in the UV, satellite measurements provide scene reflectivities that are quantified with LERs at wavelengths in the range of weak gaseous absorption. To derive LER rs from a measured radiance IM, the atmospheric path radiance Ia, transmissions T and T, and reflectance Sb for a spectral band are calculated for a molecular atmosphere, and the inversion of Eq. (5) yields

(6) r s = I s T T + S b I s ,

where Is=IM-Ia. A vast majority of scene LERs derived from satellite observations contain contributions from scattering from clouds, aerosols, or both (see Sect. 2.3 for their treatment). To characterize reflective properties of natural surfaces, many investigations have been devoted to creating global LER climatologies by selecting gridded LERs that are minimally affected by clouds or aerosols from repeated observations over a period of time (typically a calendar month). These climatologies include spectral surface LER databases constructed from the TOMS radiance measurements at 340–380 nm from 1978–1993 (Herman and Celarier1997), GOME-1 at 335–772 nm from 1995–2000 (Koelemeijer2003), SCIAMACHY at 335–1670 nm from 2002–2012 (Tilstra et al.2017), OMI at 328–499 nm from 2005–2009 (Kleipool et al.2008), and GOME-2 at 335–772 nm from 2007–2013 (Tilstra et al.2017). Intercomparisons of these spectral LERs from different satellite missions show good agreement among corresponding measurements (Tilstra et al.2017) despite differences in observation time periods, viewing illumination geometry, and footprint size. For a location on Earth, its surface is usually observed at nearly the same local solar time from a Sun-synchronous orbit, and thus the sampling of its surface BRDF is limited to a small range of SZAs. Furthermore, the selection of cloud- and aerosol-free LERs tends to favor low LER values, thus likely excluding the LERs at high VZAs. LER values of natural surfaces tend to be quite close when SZAs fall within a small range and large VZAs are excluded; hence, these LER climatologies are presented as independent of viewing illumination geometry. The low LER sensitivity to varying viewing illumination geometry (within limited ranges of SZA and VZA) indicates that natural surfaces (excluding glittering water surface) have weak anisotropy and can be treated as Lambertian surfaces. These climatological data reveal that the surface LER in the UV for snow- and ice-free areas varies within the range of 0.02–0.1 for most land and (off-glint) water surfaces, except for a few places on Earth, such as the Saharan desert and the salt flat in Bolivia, where surface LERs may exceed 0.1. These low surface LER values derived from satellite observations have been validated in field experiments (Coulson and Reynolds1971; Doda and Green1980, 1981; Feister and Grewe1995), which have found that the spectral reflectances of natural surfaces, such as the open ocean, forest, grassland, and desert, fall within the same range of satellite LER measurements. These field experiments have also demonstrated that the spectral reflectances of natural surfaces vary slowly and smoothly with changing wavelengths. The spectrally smooth GLER of natural surfaces permits accurate estimation of GLER within the UV range with measurements at two or more wavelengths, specifically the extrapolation of GLERs determined at the long (weak O3 absorption) wavelengths to estimate the GLERs at short (strong O3 absorption) wavelengths.

Based on the reflective characteristics of natural surfaces described above, the forward model for retrieval treats the reflections from a surface as Lambertian, whose reflectance is determined from the radiance measurement of the spectral band with weak gaseous absorption or is extrapolated from the weak to the strong absorption band. We use the reflection from an ocean surface as an example to illustrate the success and deficiency of the isotropic surface treatment and the GLER extrapolation, since a water surface is likely the most anisotropic surface encountered in satellite remote sensing. Figure 4a displays the GLERs of an ocean surface at the four EPIC UV bands as a function of VZA along the incident plane with the Sun at θs=15. Viewing in the specular direction (θv=15 and ϕ=180), the GLER decreases with longer wavelengths, but the reverse is true when viewing in directions 25 or greater away on either side of it. In other words, the reflection appears to be less anisotropic at shorter wavelengths. This is due to a less direct beam and thus more diffuse radiation (resulting from more photons being Rayleigh-scattered by air molecules) at the shorter wavelengths. While the reflection of a direct beam yields anisotropic outgoing radiation according to the BRDF, the diffuse radiation impinges on the surface from every possible direction of the hemisphere above, usually resulting in much less anisotropic reflected radiation, which follows the angular distribution specified by the hemispherically averaged BRDF. Figure 4b shows the spectral dependence of GLER on wavelength, illustrating that linear extrapolation of GLER at longer wavelengths (340.0 and 388.0 nm) yields highly accurate GLER estimations at shorter wavelengths (317.5 and 325.0 nm), usually with errors much less than 1 %.

The Lambertian surface treatment enables an accurate estimation of the surface radiance Is without knowledge of the actual BRDF, provided that the GLERs estimated at some (usually the weak absorbing) wavelengths can be accurately extended (linearly extrapolated) to other wavelengths. However, the paths traversed by photons reflected from a Lambertian surface differ from those from an anisotropic one, as illustrated in Fig. 5, which displays the mean optical path lengths, ms=-lnIs/τz, of EPIC band 1 as a function of altitude for two viewing illumination geometries. As shown in Fig. 5, the path lengths differ the most just above the surface, but the difference decreases with higher altitudes due to less course-altering atmospheric scattering resulting from lower air density and vanishes around 25 km above the surface. Thus, the Lambertian treatment of an anisotropically reflective surface may introduce an error, called the AMF error, in accounting for atmospheric absorption due to the difference in the photon sampling of the atmosphere. This difference is larger in the lower troposphere but becomes negligible in the stratosphere, implying that the effect of anisotropic reflection, i.e., the BRDF effect, has a larger impact on the quantification of trace gas absorption in the troposphere but a smaller one for trace gases in the stratosphere. Because the bulk O3 (∼90 %) is located in the stratosphere, the Lambertian treatment does not introduce a significant AMF error in total O3 absorption.

Figure 5Mean path lengths (ms) of EPIC band 1 reflective photons from an ocean surface (with the same BRDF described in Fig. 4) and its Lambertian equivalent surfaces. Here the mean path lengths ms, normalized by the respective geometric air mass factors (mG), are plotted as functions of altitude z for two viewing illumination geometries: one view from the direction of specular reflection at θv=15 and ϕ=180 and the other at θv=50 and ϕ=0, while the Sun is at θs=15 for both geometries.


As described above, UV reflectivities for most natural surfaces are quite low (GLER <0.1); therefore, the surface contributions Is are typically much smaller than (<10 % at 317.5 nm) the path radiance Ia (see Fig. 6). In modeling a measured radiance IM, an error in surface radiance Is is compensated for with the path radiance Ia. The uncertainty of extrapolated GLER is usually less than 1 %, corresponding to a less than 1 % error in Is and hence less than 0.1 % error in the path radiance Ia. Furthermore, the AMF error due to the Lambertian treatment of an anisotropic surface is insignificant, since the combined mean photon path lengths,

(7) m z = - ln I TOA / τ z = ( I a m a + I s m s ) / I TOA ,

contain minor contributions from surface radiance Is.

Figure 6Example path radiance, Ia, and surface radiance Is for θs=45, θv=40, and ϕ=135. Ia is the middle line in black, and Is for LER =0.1 and LER =0.94 is represented by the lower (red) and upper (blue) lines, respectively.


Natural surfaces with high UV reflectivities (GLER >0.2) are surfaces covered with snow, ice, or both. The highest GLER values are found over Antarctica and Greenland, where typical GLER values are higher than 0.9, as shown in Fig. 7. Figure 7 shows sample results of a climatological GLER database for Antarctic ice constructed from the observations of polar-orbiting instruments, including Aura OMI and SNPP OMPS, and it reveals a sizable dependence of ice GLER values on the viewing illumination geometry, indicating that the reflection from ice is significantly anisotropic. Because of the much higher surface radiance Is (e.g., Fig. 6 blue line), the Lambertian treatment of ice surface can lead to large AMF errors. However, the ice GLER varies within a small range (0.94 to 0.98), and hence ice reflection has weak anisotropy for low SZA and VZA (<70). Because the stronger O3 absorption and Rayleigh scattering at shorter wavelengths reduce the fraction of direct solar beam but increase that of the diffuse radiation reaching the surface, further weakening the BRDF effect, the error of Lambertian treatment of ice surface in the sampling of atmospheric O3 absorption is suppressed for the low SZA and VZA observations.

Figure 7Climatological Antarctic GLER values at 331 nm as functions of SZA (θs) for three viewing geometries, revealing a significant dependence of ice GLER on the viewing illumination geometry.


2.3 Particle scattering and absorption

Atmospheric particles, including clouds and aerosols, reside mostly in the troposphere and cover a large portion (∼67 % by clouds alone; King et al.2013) of the Earth's surface. Radiative transfer modeling of sunlight through a particle-laden atmosphere can be performed to quantify the TOA contributions from possible light paths, provided that the optical (scattering and absorption) properties of these particles, their amounts, and vertical distributions are specified. However, for UV remote sensing observations, the quantitative information about particles needed for radiative transfer modeling is in general not sufficiently known, precluding their explicit treatment. In this section, we describe an implicit treatment of atmospheric particles for the simulation of measured radiances with the mean photon path approximately matching that through the particle-laden atmosphere.

Atmospheric particles scatter and possibly absorb UV photons; they can thus significantly alter their paths through layers from closely above the particles down to the ground surface, usually shortening the path lengths below while lengthening those above the particles. Observing from space, the apparent effect of atmospheric particles is the enhancement of the TOA radiance contributed by backscattering from them. Since this effect is very similar to the consequence of an increased surface albedo, it is often referred to as the albedo effect. The albedo effect can be modeled by placing in a molecular atmosphere an elevated bright surface that partially covers an IFOV. This treatment is called the mixed Lambertian equivalent reflectivity (MLER) model, which is frequently employed by many algorithms for trace gas retrievals. Based on the MLER model, the TOA radiance for an IFOV is expressed as

(8) I TOA = I g ( R g , p g ) ( 1 - f c ) + I c ( R c , p c ) f c ,

the weighted sum of two independent contributions Ig and Ic. Here Ig is the radiance from the cloud-free portion of the IFOV containing a Lambertian surface of reflectivity Rg at pressure pg. Similarly, Ic is from the cloudy portion, and fc is the cloud fraction and Rc the reflectivity of the Lambertian surface at pressure pc.

The MLER model can reproduce measured radiances Im through the determination of cloud fraction fc. First, the scene LER rs at surface pressure pg is estimated using Eq. (6). If rs is less than or equal to the climatological LER value Rg (e.g., Kleipool et al.2008), this IFOV is treated as a particle-free scene (fc=0). If rs is greater than or equal to the LER value for cloud Rc=0.8 (Koelemeijer and Stammes1999; Ahmad et al.2004), this IFOV is treated as fully cloud-covered (fc=1). When rs is between Rg and Rc, the cloud fraction is inverted from Eq. (8), which yields

(9) f c = I M - I g I c - I g .

In the case of fc=0 or 1, surface LER rg or cloud LER rc is determined using Eq. (6) to ensure that modeled radiance ITOA is equal to the measurement IM. Figure 8 shows cloud fractions (fc) as a function of wavelength for several examples of particle-laden atmospheres.

Figure 8Four examples of cloud fractions (fc) derived from explicitly modeled TOA radiances for particle-laden atmospheres. The first of these is the atmosphere with a 1.5 km thick layer of C1 cloud (CLD; Deirmendjian1969) with a single-scattering albedo ω=1 and an optical thickness τ=5 at 340 nm centered at 5 km altitude (or pressure level of 545 hPa). The others are atmospheres with a 1 km thick layer of aerosols, including SLF (ω=0.996), BIO (ω=0.921), and DST (ω=0.900) aerosols (SLF, BIO, and DST models are taken from Torres et al.2007), with an optical thickness τ=1.5 at 340 nm centered at 3 km altitude (or pressure level of 703 hPa). The insets list the MLER parameters, Rg, pg, Rc, and pc, as well as the angles (θs, θv, and ϕ) that specify the viewing illumination geometry.


The radiance intensity scattered from atmospheric particles varies smoothly with wavelength without high-frequency spectral structures. For instance, the contributions to TOA radiances (ITOA) from backscattering by meteorological clouds change smoothly and slowly with wavelength (see Fig. 8, the CLD curve). The selection of Rc=0.8 facilitates close simulation of the spectral variation of clouds observed from space (Ahmad et al.2004) by the MLER model such that retrieved fc has a small spectral variation (i.e., fc nearly the same for different wavelengths) for most cloudy observations. The small and smooth change of fc with wavelengths allows its extrapolation to provide a reliable estimate of fc at shorter wavelengths from those determined at longer wavelengths.

Certain types of aerosols, such as continental aerosols containing soot, smoke from fires, mineral dust from deserts, and ash from volcanic eruptions, both scatter and absorb UV photons passing through them. Usually, aerosol absorptions cause the underlying surface (including clouds) to appear darker, more so at shorter wavelengths. The change in ITOA due to the addition of aerosols and hence the cloud fraction (fc or the surface LER, rg) are smooth in wavelength (see, e.g., Fig. 8, with smooth curves for weakly absorbing sulfate-based aerosols – SLFs, carbonaceous aerosols from biomass burning – BIO, and mineral dust – DST; Torres et al.2007). Therefore, fc (when fc>0, from Eq. 9) or rg (when fc=0, from Eq. 6) determined at longer wavelengths where atmospheric absorption is weak may be linearly extrapolated to O3-sensitive wavelengths for estimation of contributions to TOA radiance from surface reflection and particle backscattering (referred to as the rgfc extrapolation method hereafter).

The UV aerosol index (AI; Herman et al.1997; Torres et al.1998), which measures the deviation of the spectral variation of TOA radiance from that of a pure molecular atmosphere, is proportional to the spectral slope cl used in the rgfc extrapolation scheme. Algebraically, AI is calculated as the N-value (defined as −100log 10I) difference between the modeled (ITOA) and measured (IM) radiances at a wavelength λ.


Here, the modeled radiance ITOA(λ,Re) is calculated for a molecular atmosphere with an estimated reflectivity parameter Re, which may be the LER value rs or the MLER cloud fraction fc determined at a well-separated wavelength (λλ). The pair of wavelengths used for the AI calculation is in the UV spectral range with weak molecular absorption, and their separation Δλ should be sufficiently large (>10 nm) to capture the spectral contrast of Rayleigh scattering. Using IM(λ)=ITOA(λ,Rm)=ITOA(λ,Re+ΔR), since the reflectivity parameter Rm is derived from IM(λ) and ΔR=Rm-Re=clΔλ, we arrive at Eq. (11) from the definition of AI in Eq. (10). In short, the spectral slope cl is equivalent to the AI, which is significantly positive for particles (such as smoke, dust, and volcanic ash) with large absorption and slightly positive to negative for non-absorbing and weakly absorbing particles (such as clouds and sulfate aerosols). Note that for the conventional AI (also called LER AI) calculation, radiance ITOA is modeled for a Rayleigh-scattering-only atmosphere over a Lambertian surface. To capture the spectral slope of the rsfc extrapolation scheme, we switch the LER treatment with the MLER modeling of ITOA for AI calculation. The resulting MLER AI is usually higher than the corresponding LER AI when fc>0, but otherwise it can be similarly used to indicate the presence of UV-absorbing aerosol.

The MLER treatment enables the modeling of measured radiances without knowledge of the optical properties or the full vertical distributions of atmospheric particles. The accuracy of the modeled radiances at the extrapolated wavelengths depends on how closely the MLER parameter (rg or fc) follows the linear relationship among different wavelengths. In reality, the spectral dependence of natural surface reflection (rg) or particle scattering and absorption (fc) is nonlinear, though moderately as exemplified in Figs. 4b and 8. Therefore, rgfc extrapolation yields small errors in rg or fc at the extrapolated wavelengths. The radiance uncertainties associated with the rgfc extrapolation error are below 1 % for the vast majority of remote sensing observations. Higher radiance uncertainties usually occur in the presence of highly elevated or strongly absorbing aerosols. These observations may be flagged with high AI values.

In addition to the mostly small radiance errors at the extrapolated wavelengths, the MLER treatment can simulate the photon sampling of particle-laden atmospheres with a diverse range of particle types and vertical distributions. Figure 9 shows comparisons of mean photon path lengths of particle-laden atmospheres with those from the corresponding MLER treatments. These comparisons illustrate that the layer mean photon paths based on the MLER model deviate from those of the particle-laden atmospheres, mostly in the region immediately above the particles down to the underlying surface. These deviations diminish with higher altitudes where lower air density reduces the chance of photons being scattered. Since the vast majority of clouds and aerosols are in the lower troposphere (<10 km), the MLER treatment does not introduce significant AMF errors in accounting for O3 absorption, which occurs mostly in the stratosphere. This is similar to how the Lambertian treatment of surface reflection works for the estimation of total O3 absorption (see Sect. 2.2).

Figure 9Mean photon path lengths mz, normalized by the geometric AMF mG, of EPIC bands 1 and 2 as functions of altitude z for particle-laden atmospheres and their MLER treatments. See the caption of Fig. 8 for the description of aerosol characterizations, MLER treatments, and the viewing illumination geometry.


The MLER treatment relies on a few adjustable parameters, including the cloud fraction fc and cloud pressure pc, to model a vast range of conditions encountered in remote sensing of Earth's atmosphere. The cloud fraction fc, obtained directly from radiance measurements using Eq. (9), provides an estimate of the cloud amount in an IFOV. The pressure pc of the elevated Lambertian surface needs to be set at a proper level to best approximate the layer mean photon paths of a particle-laden atmosphere. As seen in Fig. 9, the optimal placement of the elevated Lambertian surface is within the particle layer, as pc located too high or too low from the optical centroid pressure (OCP; Joiner and Vasilkov2006; Vasilkov et al.2008) would make layer mean photon paths deviate further from those of the particle-laden atmosphere. The effective cloud pressures retrieved from the EPIC measurements of the O2 A band (Yang et al.2019) are usually located within the particle vertical distributions and therefore used to set the cloud pressures pc for processing EPIC observations.

The use of OCP for pc enables the MLER model to account for the measurement sensitivity change when a layer of particles is introduced into the atmosphere: enhancing the photon attenuation by absorbers inside and above the layer, while reducing them below, as the mean photon paths or AMFs from the MLER model lengthen above pc but shorten below it, as illustrated in Fig. 9. Since the MLER model captures the enhancement and shielding effects on trace gas absorption by atmospheric particles, it is widely adopted due to its simplicity for retrievals of trace gases besides O3, such as NO2 and SO2 in the troposphere. However, sizable AMF errors are prevalent for modeling tropospheric absorptions based on the MLER treatment, which usually yields significantly different mean photon paths from those of explicit treatment in the troposphere.

2.4 Inelastic molecular scattering

The scattering of sunlight with atmospheric constituents is mostly elastic; i.e., the energy and thus the wavelength of a photon remain the same before and after the interaction. But a small portion (∼4 %) of molecular scattering is inelastic, resulting in energy gain or loss of the scattered photons. Specifically, the rotational Raman scattering (RRS) from air molecules (such as nitrogen and oxygen) can alter the wavelengths of scattered photons, with UV wavelength shifts Δλ<±2 nm (Joiner et al.1995; Chance and Spurr1997; Vountas et al.1998). These inelastic scatterings cause the filling-in of telluric lines (i.e., trace gas absorption features) and solar Fraunhofer lines (also known as the Ring effect, which was first noticed by Grainger and Ring1962).

The filling-in effect is a function of wavelength and depends on the optical properties of the atmosphere, the viewing and illumination geometry, and the surface reflectivity and pressure. The filling-in effect also depends on the ISRF, especially on the instrument spectral resolution, which is the width of its ISRF, since the measured radiance of a band is a convolution of spectral radiance and the ISRF (see Eq. 1). This effect is quantified with the filling-in factors, defined as ϱ=(IRRS-IELA)/IELA, where IELA is the TOA radiance calculated assuming all molecular scattering is elastic, while IRRS includes the inelastic (RRS) contributions. To illustrate the significance of RRS, we show in Fig. 10 examples of the filling-in factors, calculated for EPIC bands using the scalar LIDORT-RRS radiative transfer code (Spurr et al.2008). Since RRS is weakly dependent on polarization, a scalar radiative transfer model, from which both IELA and IRRS are calculated without including radiation polarization, can accurately provide filling-in factors (Landgraf et al.2004; Wagner et al.2010).

The filling-in factors provide estimates of the modeling errors in ITOA when RRS contributions are neglected, and results in Fig. 10 show variations of modeling errors with different observing conditions. These errors are usually systematic for a spectral band and are between 0.5 and 1 % for measurements of EPIC bands 1 and 2. These errors are sufficiently large that corrections are required for achieving high (1 %) O3 retrieval accuracy. The filling-in factors (ϱ), modeled using a scaler code (like LIDORT-RRS), may be used to correct the results (ITOA) from vector radiative transfer codes (e.g., Dave1964; Spurr2006) that perform elastic modeling only, i.e., the RRS-corrected TOA radiance =ITOA(ϱ+1).

Figure 10Filling-in factors ϱ=IRRS-IELAIELA×100 for EPIC UV bands as a function zenith angle at two surface pressures pg=0.7 ATM and pg=1.0 ATM, with a surface albedo of rg=0.1.


3 O3 and temperature vertical profiles

As shown in Sect. 2.1, the O3 vertical distribution or profile directly affects the magnitude of a measured radiance in the spectral region with significant O3 absorption. Hence, the interpretation of radiance change due to O3 absorption requires some knowledge of its profile. In general, the retrieval of quantitative information about a gaseous absorber (such as O3 and SO2) requires a model to prescribe its vertical distribution. The skill of this model in representing the actual vertical distribution of the absorber significantly contributes to the quantification accuracy. In this section, we describe a recently developed O3 profile model for remote sensing retrieval algorithms and its improvements over the model commonly used by other total O3 algorithms.

Ozone is naturally present throughout the atmosphere, and its spatial and temporal distribution is controlled by atmospheric processes of O3 production, destruction, and transport. The O3 distribution exhibits a high abundance of O3 in the stratosphere and a minor portion (∼10 %) in the troposphere, with the peak O3 concentration occurring at a lower altitude as the latitude increases towards the poles. These characteristics are well-captured by O3 profile climatologies (e.g., Fortuin and Kelder1998; McPeters et al.2007; McPeters and Labow2012), which provide the mean and variance of O3 vertical distribution as a function of latitude and calendar month. These climatologies also reveal that the O3 profile has the highest variability in the upper troposphere and lower stratosphere (UTLS), contributing the most to the natural variations in total O3. This high O3 variability is the consequence of atmospheric movements that blend air masses with different O3 concentrations, such as uplifting of O3-poor air in the troposphere or lowering of O3-rich air in the stratosphere resulting from the rise and fall of the tropopause. Predictors of O3 profile shape, including tropopause pressure and total O3 columns, are developed to capture the dynamical influences on O3 vertical distributions, resulting in the construction of tropopause-sensitive (Wei et al.2010; Bak et al.2013; Sofieva et al.2014) and total-column-dependent (Wellemeyer et al.1997; Bhartia and Wellemeyer2002; Lamsal et al.2004; Labow et al.2015) O3 profile climatologies.

The O3 profile model for the Total Ozone Mapping Spectrometer Version 8 (TOMS-V8) total O3 algorithm combines the latitude-dependent monthly mean Labow–Logan–McPeters (LLM) climatology (McPeters et al.2007) with the latitude- and total-column-dependent annual mean climatology (Bhartia and Wellemeyer2002) to determine the O3 profile as a function of latitude, time (day of year, DOY), and total O3 column. This model has been adopted by nearly all the contemporary total O3 algorithms (e.g., Bhartia and Wellemeyer2002; Eskes et al.2005; Veefkind et al.2006; Van Roozendael et al.2006; Lerot et al.2010; Loyola et al.2011; Van Roozendael et al.2012; Lerot et al.2014; Wassmann et al.2015) owing to its capability to characterize O3 profile variation with the total column.

To improve the representation of the O3 profile, we construct both tropopause-dependent and total-column-dependent climatologies using the Modern-Era Retrospective Analysis for Research and Applications Version 2 (MERRA-2; Bosilovich et al.2015; Gelaro et al.2017) O3 record between 2005 and 2016. The total-column-dependent climatology, named M2TCO3, is more appropriate for use as the O3 profile model needed by a total O3 algorithm, as it is generally more reliable than the tropopause-dependent version in prescribing realistic O3 profiles (Yang and Liu2019).

Figure 11 compares daytime M2TCO3 (Yang and Liu2019; referred to as M2TCO3 hereafter) and TOMS-V8 profiles for 2 months and four latitude zones, illustrating the similarities and differences between the two O3 models. Both show north–south asymmetry, i.e., profiles in the Northern Hemisphere differ from those in the Southern Hemisphere for the corresponding months and latitude zones (e.g., September and 60–50 S vs. March and 50–60 N in Fig. 11), substantial seasonal variations (e.g., 60–50 S, March vs. September in Fig. 11), a strong dependence on latitude exhibiting lower altitudes of O3 concentration peaks at higher latitudes for similar total columns, and characteristic dependence on the total column, which gets smaller with a higher O3 peak altitude (e.g., March and 50–60 N in Fig. 11). Figure 11 shows good agreements of zonal mean profiles (e.g., close matches between solid black and dotted black curves in each panel of this figure) but significant differences between M2TCO3 and TOMS-V8 profiles for similar total columns. These differences are due to TOMS-V8's use of annual mean column-dependent climatology to account for profile variations with the total column throughout the year (Bhartia and Wellemeyer2002), thus ignoring the significant seasonal dependence. An additional deficiency of TOMS-V8 contributing to the differences is its inadequate representation of latitude-dependent O3 profile variation with the total column, including broad (30) latitude zones and omission of north–south asymmetry. These deficiencies are eliminated with M2TCO3, which improves the realism of O3 profile representation.

Figure 11Profile comparisons between M2TCO3 and TOMS-V8 for 2 months (March and September) and four latitude zones: 60–50 S, 30–20 S, 20–30 N, and 50–60 N. Colored solid lines represent M2TCO3 profiles, while the dotted ones are for TOMS-V8 profiles. The color of a solid line indicates the percentage occurrence of the climatological profile, and its line legend displays the mean tropopause altitude and the mean total column O3 of the profile. The solid black lines represent the downgraded M2TCO3 (i.e., the monthly zonal mean) profiles, and dotted lines are TOMS-V8 monthly zonal mean (i.e., the LLM climatological) profiles. Here pressure altitude is defined as Z*=16log10psp, where p is pressure level (in hPa) and ps=1013.25 hPa.


In short, M2TCO3 better captures the dynamical changes and spatiotemporal variations in O3 profiles with higher resolutions in total O3 column (25 DU), latitude (10), and time (monthly). Taking into account the substantial change in atmospheric O3 over a long time, M2TCO3 is more accurate in representing atmospheric O3 vertical distribution from the recent past to near future than the TOMS-V8 model, which was compiled from earlier satellite and ozonesonde data (mostly from the 1980s and 1990s; Wellemeyer et al.1997; McPeters et al.2007). Hence, we use the M2TCO3 climatology as the O3 profile model for total O3 retrieval from EPIC.

The M2TCO3 climatology contains not only mean profiles that represent the likely O3 vertical distributions, but also the modal O3 adjustment profiles that specify the probable deviations from the means. These modal profiles are determined from the O3 profile covariance statistics, as illustrated in Fig. 12, showing examples of M2TCO3 climatological O3 profiles and the associated modal profiles, which are the eigenvectors (also known as the empirical orthogonal functions or EOFs) of the profile covariance matrices. Algebraically the representation of an O3 profile X is expressed as

(12) X = X m ( v ) + k = 1 p γ k e k ( v ) ,

where Xm(v) is a climatological profile that depends on a set of variables v, which for M2TCO3 consists of the total column (Ω0), time, and location. ek(v) is the kth modal profile, γk the kth coefficient, and p the number of ek(v), with a maximum equal to the number of levels used to represent an O3 profile in the climatology. Usually, a few modal profiles are sufficient to account for the majority of profile variance. For example, in Fig. 12, the first five EOFs (panels b and e) of the covariance matrices (panels c and f) account for 80 % of profile variances (blue shaded area in panels a and d). An actual O3 profile X, which deviates invariably from the mean Xm, can be accurately represented using Eq. (12) with a small number of expansion coefficients γk. Much like the mean the profile Xm represents the most probable vertical distribution of O3: the modal profiles, {ek,k=1}, describe the most, the second most, and so on likely vertical patterns of deviations from the mean profile. Each modal profile describes a rearrangement, like shifting, shrinking, or broadening, of the mean profile without substantially changing the total column. With these modal profiles constraining how a profile can be adjusted, the retrieval algorithm can exploit the O3 profile information contained in multi-spectral measurements to improve the O3 profile representation by determining one or more linear expansion coefficients {γk,k=1}. Note that for most total ozone algorithms, the O3 profile representation is limited to the climatological mean only, equivalent to restricting γk=0 for all k in Eq. (12).

Figure 12Examples of M2TCO3 climatological profiles for the southern midlatitude zone in March (a) and the northern midlatitude zone in September (b), the associated correlation matrices (c, f), and the corresponding modal O3 profiles (b, e). The blue shaded areas in panels (a) and (d) are within 1 standard deviation of the mean. The correlation matrices in (c) and (f) are standardized (i.e., diagonal element normalized to 1) covariance matrices. The five modal profiles in (b) and (e) are the first five ordered eigenvectors (also known as empirical orthogonal functions or EOFs) of the corresponding covariance matrices, with percentages of the profile variance explained by the EOFs displayed in the line legends. The text box in each panel displays the average tropopause altitude (in km) and the average total O3 column (in DU) for the climatological profile.


The total column is a good predictor of an O3 profile that is especially accurate for the shape in the stratosphere but less so in the troposphere. Tropospheric O3 exhibits characteristic spatiotemporal distribution, which is captured in the MERRA-2 tropospheric O3 climatology (Yang and Liu2019). To better represent the O3 profile, the tropospheric part of a column-dependent M2TCO3 profile, Xm, is scaled with the ratio of the MERRA-2 climatological tropospheric column to the tropospheric column integrated from the downgraded M2TCO3 profile (see Fig. 11 for sample M2TCO3 and downgraded M2TCO3 profiles). In other words, the profile Xm in Eq. (12) has its tropospheric part tied to the spatiotemporally varying climatological tropospheric column, to which the tropospheric column of the mean Xm profile (obtained by averaging over the different column amounts) is matched.

In addition to knowledge of profiles of light-absorbing trace gases, such as O3 and SO2, radiative transfer modeling of measured radiance requires knowledge of the atmospheric temperature profile because the absorption cross-sections of these trace gases significantly depend on temperature. For total O3 retrieval from EPIC, this knowledge is taken from the temperature profile climatology created from MERRA-2 data together with the ozone profile climatology (Yang and Liu2019). This temperature climatology provides mean temperature profiles corresponding to the climatological O3 profiles, capturing the dependence of the temperature profile on season and location, as well as the variation of temperature with the O3 profile. It is an improvement over the TOMS-V8 temperature profile climatology, which provides latitude- and month-dependent temperature profiles, but without accounting for the strong correlation between temperature and O3 profiles.

4 Inversion technique

Section 2 describes algorithm physics treatments of interactions of solar radiation with atmospheric particles and surfaces to enable RT modeling of photons traversing through a molecular atmosphere to reproduce the measured TOA radiances with photons that follow the paths similar to those through the actual atmosphere and therefore establish the relationship between spectral measurements and the atmospheric state, as well as surface reflectivity and instrumental parameters. At its core, the RT modeling sets up a forward mapping from the vertical distributions of gaseous absorbers and the surface reflectivity parameters to measured TOA radiances. The retrieval of gas absorbers, such as O3 and SO2, is the inverse of this mapping, i.e., to find their vertical distributions and the surface reflectivity parameters for which forward modeling closely reproduces the measured TOA radiances. However, this inverse mapping is inherently an ill-posed problem, as the solution is not unique; i.e., more than one set of profiles and surface parameters can yield the same measurements. This problem is made worse with measurement uncertainties, which expand the profile and surface combinations that can reproduce, within error bars, the measured spectra.

For successful inversion, analytical constraints are placed on the profiles of gas absorbers as well as the spectral variations of ground reflectivity and atmospheric particle (aerosol and cloud) backscattering. For O3 retrieval, Eq. (12) embodies the profile constraint, while the MLER model with rsfc extrapolation regulates the surface reflection and particle backscattering. These constraints control the dimension of the inverse mapping space and manifest themselves as the retrieval (i.e., adjustable) parameters, which, in the case of O3 retrieval, consist of total O3 column Ω, a number (p) of modal expansion coefficients {γk,k=1p}, surface LER (rs) or cloud fraction (fc), and a number (q) of polynomial coefficients {cl,l=1q} of the rsfc extrapolation. The set of adjustable parameters forms the state vector (x) whose length (n) is the dimension of inverse mapping space. Proper selection of adjustable parameters by limiting the number of modal coefficients (p≥0) and polynomial coefficients (q≤1) ensures that the inverse problem is well-posed and simultaneously maximizes the amount of information collected from the spectral measurements. Here p=0 indicates no modal expansion, equivalent to restricting the profile to a climatological column-dependent O3 profile, and q=0 for the spectral invariant reflectivity parameter.

4.1 Exact solution

Conceptually, the inversion is to find the state vector (x) that satisfies a set of m simultaneous equations, {Δyi=0,i=1m}, one for each spectral band difference, Δyi=lnIM(λi)-lnITOA(x,λi), between the radiance measurement IM and the forward modeling ITOA: here λi the wavelength that characterizes the ith (1im) spectral band and Δyi the residual of this band. In matrix form, the m simultaneous equations can be expressed as

(13) Δ y = 0 ,

where Δy is residual column vector {Δyi,i=1m}. Since the forward mapping ITOA(x) is a nonlinear function of the state vector x and has no analytical inverse, the solution to Eq. (13) is usually obtained iteratively. The iteration is started with an initial (i.e., iteration number L =0) state vector xL to linearize the equation between residuals and the state vector.

(14) Δ y i = ln I M λ i - ln I TOA x L , λ i - j = 1 n x j - x L j ln I TOA x , λ i x j | x = x L ,

where xj and xLj are the jth components of x and xL, respectively, Δxj=xj-xLj the jth components of state adjustment vector, and Kij=lnITOAx,λixj|x=xL the Jacobian, also known as the weighting function for the retrieval parameter xj at spectral band λi. The m residual elements, each written in Eq. (14), can be expressed in matrix form as

(15) Δ y = Δ y L - K Δ x ,

where ΔyL is the column vector {lnIM(λi)-lnITOA(xL,λi),i=1m}, Δx=x-xL the state adjustment vector, and K the m×n Jacobian matrix with the {Kij,i=1m,j=1n} as its elements. Putting Eq. (15) into Eq. (13) yields

(16) Δ y L = K Δ x ,

which may be solved exactly (under strict conditions) to determine the state adjustment vector Δx. After each iteration, the linearization state vector is updated to

(17) x L + 1 = x L + Δ x .

The final state x is found when the iteration converges, i.e., when the absolute change in state vector Δx is below a threshold.

The linear equation (Eq. 16) may be solved exactly only when the number of measurements is equal to the number of retrieval parameters (i.e., m=n) and the Jacobian matrix K is invertible (i.e., non-singular matrix), as exemplified in the well-known TOMS-V8 total O3 algorithm (Bhartia and Wellemeyer2002). The TOMS-V8 algorithm determines the two-component state vector, x={Ω,rsorfc}, from radiance measurements of two spectral bands: one with low O3 sensitivity to estimate the MLER parameter (rs or fc) and the other with high O3 sensitivity to derive total O3 column Ω. However, few other algorithms adopt this inversion method, since it requires m=n and K to be a non-singular matrix. Even if both these conditions are met, inverting Eq. (16) to obtain exact solutions tends to enhance the impact of measurement uncertainties (noises) on the retrieved results, as in cases that K matrices are nearly but not quite singular. These cases occur when the spectral variation of a Jacobian has some similarity or a high degree of correlation with that of another retrieval parameter, leading to algorithm difficulty in distinguishing two retrieval parameters corresponding to the two Jacobians, thus yielding unstable retrieval results, such as in the case of simultaneous retrieval of total O3 and SO2 columns from EPIC UV measurements.

4.2 Direct fitting

Since spectral measurements have errors and mn in general, the inversion is achieved by finding the solution x that minimizes the cost function


where Sϵ is the measurement error covariance matrix, with its ith diagonal element equal to μi2. Here μi is the fractional standard deviation of the radiance error of the ith band. In the case of independent measurement error, i.e., no error correlation between different spectral bands, Eq. (18) can then be explicitly written as Eq. (19), which is the formulation of the least-squares method.

The minimization of the cost function Υ can be started by linearizing the residuals with an initial (i.e., iteration number L =0) state vector xL. Substituting Δy (Eq. 15) into Eq. (18), we minimize this cost function to obtain the state adjustment vector

(20) Δ x = ( K T S ϵ - 1 K ) - 1 K T S ϵ - 1 Δ y L = G DF Δ y L ,

which is the solution of linear weighted least-square regression. Here, GDF=(KTSϵ-1K)-1KTSϵ-1 is the direct fitting (DF) gain matrix.

This procedure of iterative minimization of the difference between measurements and modelings to determine the bulk parameters is called the direct vertical column fitting (DVCF) algorithm. The DVCF algorithm is quite general and valid for both discrete wavelength and hyperspectral measurements, as well as for different types of retrieval parameters, such as MLER parameters, layer partial columns of various absorbing trace gases, and their total vertical columns. This algorithm has been applied to retrievals of total O3 vertical column (Joiner and Bhartia1997; Yang et al.2004; Lerot et al.2014), the combination of total O3 and SO2 vertical columns (Yang et al.2007, 2009a, 2013), the combination of O3 and altitude-resolved SO2 vertical columns (Yang et al.2009b, 2010), and stratospheric and tropospheric NO2 vertical columns (Yang et al.2014). This algorithm is named DVCF to contrast with the DOAS (differential optical absorption spectroscopy) method (Platt2017), which derives a slant column and then uses an air mass factor (AMF) at a single wavelength (λ0) to convert it to a vertical column.

In general, the DVCF algorithm works well when the changes in radiance measurements responding to changes in the state vectors are significantly different between any two retrieval parameters, i.e., that columns of K, which are the Jacobians of a retrieval parameter at different wavelengths, exhibit significantly different spectral dependence from one another. This is usually true for any two bulk retrieval parameters over a sufficiently broad spectral range, such as total O3 column (Ω) and an expansion coefficient (γk) of differential profile ek (see Eq. 12), or the SO2 vertical column and its layer altitude. With measurements from a broad spectral range, the DVCF algorithm can discriminate subtle spectral features contained in hyperspectral measurements to enhance the retrieval accuracy (e.g., Yang et al.2009b, 2010). Besides contrasting with the DOAS method, the name DVCF emphasizes the vertical column because this algorithm is usually not suitable for traditional profile retrieval due to the high similarity of partial column Jacobians between adjacent layers and hence the difficulty in distinguishing their partial columns.

4.3 Optimal estimation

In many cases, such as sparse spectral sampling or a narrow spectral range, the performance of the direct fitting inversion method may decline as a result of limited information contained in the spectral measurements. For stabilizing the retrieved results, the inversion process can be regulated with an additional constraint, which is frequently based on a priori knowledge of the retrieval parameters. Algebraically, adding an a priori constraint to Eq. (18) yields a new cost function

(21) Υ ( x ) = Δ y T S ϵ - 1 Δ y + ( x - x a ) T S a - 1 ( x - x a ) ,

where xa is the a priori state vector and Sa the a priori state vector covariance matrix. The first term on the right-hand side (rhs) of Eq. (21) strives to diminish the difference between measured and modeled radiances, performing the same function as the direct fitting retrieval, while the second rhs term seeks to reduce the deviation of retrieved x from the a priori xa. This a priori constraint effectively stabilizes the retrieval by guiding the state vector adjustment when the measurements contain little information to differentiate the contributions from different components of the state vector. Using the optimal estimation (OE) technique (Rodgers2000) to minimize the cost function, Eq. (21) yields the a posterior state adjustment vector at the Lth iteration:

(22) Δ x = S a - 1 + K T S ϵ - 1 K - 1 K T S ϵ - 1 Δ y L + S a - 1 Δ x aL

(23) = S a - 1 + K T S ϵ - 1 K - 1 K T S ϵ - 1 KG DF Δ y L + S a - 1 Δ x aL

(24) = Δ x aL + S a K T KS a K T + S ϵ - 1 Δ y L - K Δ x aL ,

where ΔxaL=xa-xL. Inserting In=(KTSϵ-1K)(KTSϵ-1K)-1, an n×n identity matrix, in front of the term KTSϵ-1ΔyL in Eq. (22) yields Eq. (23). At iteration L =0, a state vector close to the actual one is sought to be the initial state vector x0, and a frequent selection is the a priori state vector: x0=xa. This is a more robust inversion scheme that works for m>n, m=n, and m<n. Equation (24) is often used in the case of m<n, as the inversion deals with an m×m (i.e., a smaller) matrix.

Equation (23) describes the difference between the current and previous state vectors Δx as a combination of the direct fitting solution GDFΔyL (see Eq. 20, which is derived without any a priori constraint) and the difference between the a priori and the previous state vectors ΔxaL weighted by matrices KTSϵ-1K and Sa-1, respectively. For the state vector component with a strong a priori constraint, i.e., a small variance in Sa, the retrieved result gravitates towards the value of the a priori state vector, while for the one with a weak constraint, i.e., a high variance in Sa, its retrieved value is primarily determined from the measurements.

The variance of a retrieved parameter is equal to the corresponding diagonal element of the covariance matrix (Sa-1+KTSϵ-1K)-1 (see Eq. 23) and is thus less than or equal to the corresponding a priori variance in the a priori Sa matrix. In other words, the change magnitude of a retrieval parameter at each iteration is usually smaller than its a priori standard deviation. Consequently the OE method can be used as an inversion scheme to ensure retrieval stability and preserve the dependence of the retrieved results on the measurements through a careful construction of the a priori covariance matrix Sa. To further reduce the dependence on the a priori state vector, it is updated at each iteration with the linearization point, setting xa=xL, and hence Eq. (22) becomes

(25) Δ x = ( S a - 1 + K T S ϵ - 1 K ) - 1 K T S ϵ - 1 Δ y L = G Δ y L ,

where G=Sa-1+KTSϵ-1K-1KTSϵ-1 is the optimal estimation gain matrix. This setting floats the anchor point of the retrieval, allowing the measurements to drive the iteration to its final state, with the a priori covariance to limit the deviation from the anchor.

By relaxing the a priori constraints through increasing the diagonal terms (i.e., the variances) of Sa such that Sa-10, G becomes GDF and Eq. (22), as well as Eq. (25), becomes Eq. (20). In other words, the direct fitting inversion is a special case of the OE inversion scheme, which is more appropriately called the regulated direct fitting inversion. Using the knowledge of their variances (Sa) to limit some of their ranges while allowing others to change freely, the DVCF algorithm with the regulated inversion scheme is suitable for retrieving multiple parameters from discrete measurements . It is applied to EPIC UV observations for simultaneous O3 and SO2 retrievals.

5 Retrieval from EPIC UV bands

EPIC has four UV channels (see Fig. 2), referred to as B1, B2, B3, and B4, characterized by wavelengths λ1=317.5 nm, λ2=325.0 nm, λ3=340.0 nm, and λ4=388.0 nm, respectively. The radiance measurements from shorter UV channels, EPIC B1 and B2, are sensitive to both O3 and SO2 absorptions (see Fig. 13), containing information that allows the retrieval of total O3 and SO2 vertical columns, provided that the reflectivity of the underlying surface is known. This knowledge is obtained from the radiance measurements of EPIC B3 and B4, the longer wavelength channels. These channels provide information about the surface reflection and particle backscattering and have very low sensitivities to O3 and SO2 absorption such that changes in O3 and SO2 amounts result in little difference in the radiance measurements of these two bands. The reflectivity determined from B3 and B4 is used to estimate the reflectivity at the shorter wavelength (O3 sensitive) channels, accomplished with the rsfc extrapolation scheme (see Sect. 2.3). The reflectivity spectral slope cl of this extrapolation is proportional to the AI (see Eq. 11). The reflectivity parameter (R) is either the LER value rs estimated from Eq. (6) or the MLER cloud fraction fc from Eq. (9) depending on the value of fc (R=rs when fc=0, R=fc when fc>0), and its spectral slope is calculated as cl=(R4-R3)/(λ4-λ3).

Figure 13EPIC bandpass-averaged cross-sections σ for O3 and SO2 at 280 K and their ratio, ρ=σSO2/σO3.


In this section, we describe the application of the DVCF algorithm to EPIC UV measurements, the scheme to solve the difficulty arising from the non-coincidence among the different EPIC spectral observations, and examples to illustrate the success of this application.

5.1 Reflectivity correction by spatial optimal estimation (SOE)

The estimation of the O3 column from EPIC radiance measurements requires accurate reflectivity information on the underlying surface, which is extrapolated from the reflectivity determined at the longer wavelength bands (B3 and B4), but the uncertainty of this extrapolation becomes large due to EPIC's asynchronous spectral measurements. Unlike most spaceborne UV instruments which provide coincident measurements from different spectral bands, EPIC takes the spectral images sequentially, separated by a time delay of ∼30 s between adjacent UV bands. Due to the Earth's self-rotation and spacecraft jitter, different spectral images record slightly different (i.e., rotated) sunlit hemispheres. The geolocation procedure of EPIC (Blank2019) aligns different spectral images and further refines the band-to-band registration using the image correlation technique (Yang et al.2000), which is estimated to provide better than 0.1 pixel (a pixel refers to an IFOV) registration accuracy for EPIC bands. Despite this high registration accuracy, rsfc extrapolation (see Sect. 2.3) becomes less accurate for a significant fraction of EPIC IFOVs as substantial reflectivity changes may occur with small shifts in viewing and solar zenith angles since near the direct backscattering direction the particle scattering phase functions have a high angular sensitivity and the shadow areas of structured scenes change nonlinearly with viewing illumination geometry. This difficulty is unlikely to improve even with better alignment and requires a new approach to correct the extrapolated reflectivity.

The basic idea to obtain a more accurate reflectivity at an O3-sensitive band is to derive it from the radiance measurement of this band with an optimally estimated total O3 column from the nearby O3 distribution. This O3 estimation is attainable because an actual spatial distribution of the total O3 column is a smooth function of geolocation and exhibits a high degree of close-range correlation (Liu et al.2009). Algebraically, the spatial optimal estimation (SOE) method finds the reflectivity (RB) at EPIC band B by minimizing the cost function that embodies the a priori knowledge of RB and O3 spatial distribution. The first part of cost function supports a smooth (i.e., homogeneous) O3 distribution, while the second part penalizes the difference between RB and its a priori value, which is the extrapolated reflectivity (RE) from the longer wavelength EPIC bands. Hence the cost function is written as


subject to the measurement constraints {IM(λB,i)=ITOA(Ω(i),RB(i),λB), and IOFV index i=1…n, the size of the IFOV group}, which is linearized to become

(28) Ω = Ω R E + ( ln I TOA R B | R B = R E / ln I TOA Ω | Ω = Ω ( R E ) ) R B - R E Ω R E + Ω R B | R B = R E R B - R E = Ω R E + S R B - R E ,

where S=ΩRB|RB=RE. The IFOV index i is dropped in Eq. (28) without losing clarity. Here j is also an index labeling the pairing (or other) IFOV in the group, and wt(i,j) is the weighting factor that depends on the distance between the ij pair. 〈Ω〉 is the average O3 column for the group. Given RE, which is the band-B reflectivity extrapolated from the longer wavelength bands, the total O3 column Ω(RE) is retrieved from band-B radiance measurement using the exact solution method (see Sect. 4.1), and the associated O3 profile is the column-dependent M2TCO3 climatological profile Xm(Ω). The equation of measurement constraint (Eq. 28) describes a positive (since S>0 usually) linear relationship between total O3 column Ω and the surface reflectivity (in the neighborhood of RE), increasing RB requires more O3 absorption to maintain IM=ITOA.

Minimizing only the first rhs term of Eq. (26) leads to the same O3 column for all the IFOVs (i.e., {Ω(i)=Ω,i=1n}), while minimizing only the second term makes RB=RE for each IFOV. The constants α and β are weights to respectively emphasize the smoothness of O3 spatial distribution and the closeness of reflectivity between extrapolation and estimation. In the SOE scheme, weights are α=0 and β=1 for the traditional O3 retrieval, also referred to as independent-pixel retrieval, while for optimized retrieval, equal weights α=β=0.5 are used.

For optimized retrieval, the minimization of the cost function Υ (Eq. 26) can be accomplished by iteratively finding RB(i) to minimize each component Υi. The solution RB(i) that minimizes Υi is found by solving this equation:

(29) Υ i R B ( i ) = β R B ( i ) - R E ( i ) R E 2 ( i ) + α j = 1 n wt ( i , j ) ( Ω ( i ) - Ω ( j ) ) S i Ω 2 = 0 ,

which yields

(30) R B ( i ) = R E ( i ) - α R E 2 ( i ) S i ( n Ω i , R E ( i ) - j = 1 n Ω ( j ) ) β Ω 2 + α n R E 2 ( i ) S i 2 .

From Eqs. (29) to (30), only the n nearby IFOVs are included, i.e., wt(i,j)=1 for ij separation within a few (<4) adjacent IFOVs; otherwise, wt(i,j)=0. At the start of iteration, {Ω(j)=Ω(j,RE),1n}, and they are then updated using Eq. (28) with RB(i) from Eq. (30) for the next iteration, which stops until changes in {RB(i),i=1n} become sufficiently small. In practice, no more than a couple of iterations are needed to reach convergence.

Figure 14 shows an example of simultaneous retrieval from the IFOVs of an EPIC hemispheric view using the SOE method. The high-variability O3 map (Fig. 14e) from the independent-pixel retrieval contains many artifacts (high spikes and low dips in O3 columns), which are substantially reduced using the SOE method, resulting in a much more realistic (smooth) O3 map (Fig. 14f). The O3 differences (ΔΩ) between optimized and independent-pixel retrievals (see Fig. 14c, d, and h) illustrate the quantitative improvements, with a small mean O3 difference (mean ΔΩ within ±0.5 DU) and a sizable reduction in O3 noise level (standard deviation of ΔΩ≈7 DU). The corresponding reflectivity corrections are quite significant at ∼0.02 on average, with a maximum of ∼0.1 deviation from the rsfc extrapolations.

Figure 14Retrieved O3 from EPIC measurements of bands B1, B3, and B4 on 3 December 2015. (a) Optimized (i.e., α=0.5,β=0.5) O3 map based on the SOE method; (b) a comparison of optimized (orange) and independent-pixel (blue, α=0,β=1) O3 along the horizontal line (left to right) across the middle of the O3 map in (a); (c) the O3 difference map: ΔO3=O3(optimized)-O3(independentpixel); (d) the O3 difference along the horizontal line across the middle of the map in (c); (e) a zoom-in of the independent-pixel O3 map; (f) the optimized O3 corresponding to the rectangle in (a); (g) cloud fraction fc corresponding to (e) and (f); (h) O3 difference (optimized independent pixel), with a zoom-in corresponding to the rectangle in (c).

In summary, the SOE method performs single-band (B1 or B2) multiple-IFOV (or image-based) retrieval, yielding reflectivity (R) and total O3 column (Ω), with the associated profile determined by the O3 model (Eq. 12) that retains only the column-dependent M2TCO3 climatological profile Xm. The a priori knowledge of the O3 distribution, which is spatially smooth, provides the extra information to correct the initial reflectivity estimation extrapolated from the characterization based on the longer wavelength bands.

5.2 Total O3 column

Radiance measurements of EPIC B1 and B2 radiances have O3 profile sensitivity, which is higher at B1 than at B2, especially at high zenith (SZA, VZA, or both) angles (as illustrated in Fig. 3). Compared to the measurement of a single O3-sensitive band, both bands jointly provide more information that allows the refinement of climatological representation of the O3 profile. This refinement is performed by adjusting the climatological profile with the most probable modal profile (e1, see Eq. 12) so that both B1 and B2 yield the same total O3 column.

For retrieval from EPIC, the full state vector to be inverted is x={Ω0,γ1,Ξ,R1,R2}, where Ω0 is the total O3 column, γ1 the O3 profile adjustment factor, Ξ the total vertical SO2 column, and R1 and R2 the MLER parameters at EPIC B1 and B2. The regulated direct fitting of EPIC B1 and B2 radiances is applied to obtain retrieved full state vector x .

For each IFOV of EPIC, the O3 vertical column is estimated first assuming there is no SO2. The iteration starts with an initial state vector x0={Ω0=Ωc,γ1=0,Ξ=0,R1=R1S,R2=R2S}, where Ωc is the climatological total column selected from the M2TCO3 climatology based on time and location. R1S and R2S are the corrected MLER parameters at B1 and B2, respectively, using the SOE method (see Sect. 5.1).

Since EPIC radiance measurement errors between any two bands are not correlated, the measurement error covariance matrix is diagonal: Sϵ=diag(σB12=0.003452,σB22=0.003452), as estimated from the random errors of the radiance (IM) measurements (see Sect. 6.2).

There is no correlation among retrieval parameters: total O3 column (Ω0), the deviation (ω1) of the O3 profile from the mean, SO2 column (Ξ), and the MLER parameters R, except between R1 and R2. The diagonal elements of the a priori covariance matrix are Sa=diag(εΩ02=102 DU2, εγ12=22 DU2, εΞ2=0.00012 DU2, εR12=0.0012, εR22=0.0012). The off-diagonal elements are equal to zero, {Sa(i,j)=0, when ij}, except for the elements associated with R1 and R2, which may be set at Sa(4,5)=Sa(5,4)=0.98εR1εR2=0.00992, representing a high degree of correlation (0.99) between R1 and R2. This Sa essentially limits the adjustments at each iteration: |ΔΩ0|εΩ0(10 DU), |Δω1|εγ1(2 DU), |ΔΞ|εΞ(10-4 DU), |ΔR1|εR1(0.001), and |ΔR2|εR2(0.001). The strong constraint on SO2 ensures that its column Ξ never deviates far (>0.01 DU) from its initial value of 0 during the iteration, essentially enforcing an SO2-free retrieval. The strong constraints on R1 and R2 also ensure that they remain nearly the same as their initial values R1S and R2S. The constraints on O3 parameters are quite loose. Especially towards the convergence of the iteration, the absolute adjustment of each component is much smaller than the corresponding standard deviation, i.e., the square root of the corresponding diagonal element of Sa.

With the setup of error and a priori covariance matrices Sϵ and Sa, the initial state vector x0 is updated (Eq. 17) iteratively using Δx from Eq. 25 until the exit of the iteration when |ΔΩ0|<0.5 DU and |Δγ1|<0.5 DU. The retrieved total O3 column (Ω) is obtained by integrating the profile X=Xm(Ω0)+γ1e1(Ω0). In processing EPIC data, the initial O3 column Ωc of x0 for an IFOV may be set to the column Ω of a previous (or nearby) IFOV to improve the speed of convergence of the iteration.

Here we list the algorithmic procedure (see the flowchart in Algorithm 1), titled EPIC total ozone retrieval, that is applied to each EPIC level-1b (L1B) granule, which contains spectral measurements, as well as geolocation and angular information for all the IFOVs of a snapshot of the sunlit side of the Earth, to produce the level-2 (L2) O3SO2AI product. The contents of a sample O3SO2AI granule are displayed in Fig. 15, including total O3, LER, and AI, which are respectively shown in Fig. 15c, e, and f. For comparison, the MERRA-2 assimilated O3 total columns, interpolated to the time and location of EPIC IFOVs, are included in Fig. 15d, their differences O3(EPIC)−O3(MERRA-2) in Fig. 15g, and the histogram in Fig. 15h. This comparison reveals excellent agreement between MERRA-2 and EPIC total O3, showing nearly identical O3 spatial distributions with similar highs and lows. Quantitatively, the differences for samples with VZA 70 are characterized by a low mean offset (μ(EPIC)=-0.2 %) and a narrow spread (standard deviation σ(EPIC)=2.52 %). Figure 15 includes the intermediate results of the EPIC total O3 processing (see the procedure in Algorithm 1), showing the total O3 columns retrieved from B1 in panel (a) and from B2 in panel (b) using the SOE method. Both B1 and B2 total columns closely resemble the MERRA-2 (panel d) and EPIC (panel c) total O3 fields, with difference statistics showing slightly worse offsets (μ(EPIC B1)=0.25 % and μ(EPICB2)=-0.41 %) and higher standard deviations (σ(EPIC B1)=2.68 % and σ(EPIC B2)=2.68 %). The improved agreement with MERRA-2 is significant, reducing the B1 O3 spread by σ2(EPICB1)-σ2(EPIC)=0.9 % (or 2.8 DU) and the B2 O3 spread by a similar amount. These better agreements are consistent over time and location, substantiating the improved retrieval with both O3-sensitive bands over a single one, which is adopted by the TOMS-V8 algorithm. Since the MERRA-2 O3 field from the assimilation of independent measurements of the Aura OMI and Aura Microwave Limb Sounder (MLS) provides highly realistic spatiotemporal O3 representation, the smaller spread between the two-band (B1 and B2) EPIC and MERRA-2 total O3 columns indicates that the inclusion of more O3-sensitive bands enables more accurate retrievals.

Figure 15An L2 O3SO2AI granule contains the total O3 vertical columns (c), LER at 340 nm (e), and AI (f) retrieved from EPIC UV measurements at 03:53:57 UTC on 3 April 2017. (a) Total O3 column (referred to as B1 total O3 column) retrieved from EPIC B1, B3, and B4. (b) Total O3 column (referred to as B2 total O3 column) retrieved from EPIC B2, B3, and B4. (c) Total O3 from all four bands. (d) Coincident MERRA-2 total O3 columns. (g) The total O3 difference: O3(EPIC)−O3(MERRA-2). (h) The histogram of the O3 differences with SZA 70, i.e., samples within the circle in (g), with a mean difference μ(EPIC)=-0.20 % (or −0.35 DU) and a standard deviation σ(EPIC)=2.52 % (or 7.4 DU). Similarly, the O3 difference, O3(EPIC B1)−O3(MERRA-2), has a mean of μ(EPIC B1)=0.25 % (or 1.03 DU) and a standard deviation σ(EPIC B1)=2.68 % (or 7.9 DU), and O3(EPIC B2)−O3(MERRA-2) has a mean μ(EPICB2)=-0.41 % (or −1.08 DU) and a standard deviation σ(EPIC B2)=2.68 % (or 7.8 DU).

5.3 Volcanic SO2

EPIC B1 and B2 radiances respond to both O3 and SO2 absorptions, but with very different (see Fig. 13) sensitivities: SO2 is more than twice as UV-absorbent as O3 at B1; in contrast, it is significantly less at B2 and about 70 % as absorbent as O3. Consequently, the estimate of O3 absorption signals at these two bands would result in an error due to the presence of SO2 in the atmosphere: 1 DU of SO2 would usually yield more than 2 DU O3 error at B1, but only about 0.7 DU error at B2. This big difference in absorption sensitivities facilitates the detection of SO2 in the atmosphere. Given a radiance signal-to-noise ratio (SNR) of 290:1, the theoretical minimum detectable level of SO2 enhancement is ∼0.5 DU in the upper troposphere and above. However, it is difficult to distinguish SO2 at this minimum level from other changes, such as the O3 profile or surface spectral reflectance, since they can induce similar changes in the measured radiances. This difficulty is increased by EPIC's asynchronous spectral measurements, which may yield spectral variation similar to the response to adding SO2 in the atmosphere. Consequently, low levels of SO2 elevation cannot be reliably detected in EPIC observations. For significant SO2 elevations, typically those from volcanic eruptions, B1 O3 is much higher than B2 O3 (i.e., Ω12) from the total O3 retrieval (described in Sect. 5.2). Adjusting the O3 profile shape or changing the spectral reflectance of the underlying surface usually cannot eliminate this large O3 discrepancy between the two bands. Therefore, a high positive value of ΔΩ can be used to flag the presence of SO2. Furthermore, a volcanic plume usually occupies a contiguous area with a limited spatial extent. Thus, ΔΩ and Ω1 enhancements resulting from volcanic SO2 plume occur over a large group of connected IFOVs instead of isolated or a small group of disconnected IFOVs. Based on these characteristics of volcanic SO2 plumes, we next describe an algorithmic procedure to flag IFOVs with SO2 enhancements.

For reliable SO2 detection, the following procedure is applied to identify the presence of SO2 in an IFOV. First, IFOVs of likely SO2 elevations are flagged through spatial analysis of the differential O3 field (i.e., ΔΩ=Ω1-Ω2, EPIC B1 and B2 O3 difference), accomplished through contour mapping to find closed areas of local ΔΩ enhancements, i.e., areas within closed contours with Ω1 considerably higher (≥7 DU) than the Ω2 values (see, e.g., Fig. 16b). The IFOVs within this ΔΩ contour likely have SO2 elevation around 5 DU or above. Next, contour mapping of Ω1 is performed to find the longest closed contour line in the area that extends 150 pixels off the extrema of the ΔΩ contour (i.e., an image rectangle with a minimum of 300×300 IFOVs or 3000×3000 km2 that covers the ΔΩ contour). Within this closed Ω1 contour, IFOVs with likely SO2 enhancements are flagged when Ω12. This flagging is then extended to the adjacent areas outside the two contours to identify IFOVs with possible SO2 contamination. For most volcanic plumes, these two contours overlap each other greatly. Including areas within the Ω1 contour and the adjacent outside regions is designed to capture plumes with lower SO2 elevations.

Figure 16EPIC observation of the volcanic plume on 23 June 2019 from the previous day's eruption of Raikoke volcano (represented by Δ in each panel) in the central Kuril Islands of Russia. (a) B1 O3 column (Ω1) from EPIC total ozone retrieval and the elevated O3 contour. (b) B1 and B2 O3 column difference (ΔΩ=Ω1-Ω2) and elevated ΔΩ contour. (c) Vertical O3 column from EPIC total SO2 retrieval (see Algorithm 2). (d) Vertical SO2 column from EPIC total SO2 retrieval. (e) SO2 vertical column retrieved from a series of eight consecutive EPIC observations of the Raikoke plume, represented by a 1.5 km thick GDF layer centered at an altitude of 13 km above sea level.

Once detected, the SO2 quantification follows the DVCF retrieval with the initial state and a priori covariance setting described next. For the IFOV identified with SO2 contamination, the initial O3 values are spatially interpolated from the background Ω1 field, γ1=0, and initial SO2 column Ξ0=Ω1-Ω2, integrated from a vertical profile specified by a generalized distribution function (GDF; Yang et al.2010) with a width and a center altitude appropriate for the plume. The corresponding elements of the a priori covariance matrix are Sa=diag(εΩ02=102 DU2, εγ12=22 DU2, εΞ2=Ξ02 DU2), i.e., the variances associated with O3 are the same as those for total O3 retrieval, while the SO2 column variance is equal to the square of the initial SO2 estimate (Ξ0), which is a weak constraint to allow Ξ to change freely responding to the measurement. Other retrieval settings are kept the same as in the total O3 retrieval described in the previous section. We list the complete algorithmic procedure in Algorithm 2, titled EPIC total SO2 retrieval, for SO2 detection and quantification from EPIC UV observations.

The algorithmic procedure (listed in Algorithm 2) applies to regions where EPIC observes volcanic plumes to produce the EPIC volcanic SO2 product. Figure 16 illustrates the detection and quantification of volcanic SO2 from EPIC observations. Spatial analyses (i.e., contour mappings) of the intermediate results (Fig. 16a and b) of the total O3 processing (see Algorithm 1) provide reliable detection of SO2 elevations. The SO2-flagged IFOVs are then processed with the DVCF algorithm to retrieve total vertical O3 and SO2 columns simultaneously, with results shown in Fig. 16c and d, respectively. Comparison of the two O3 fields in Fig. 16 shows that the initial O3 elevations (Fig. 16a) due to the presence of SO2 are nearly entirely removed in the final O3 field (Fig. 16c), demonstrating that the combined retrieval of O3 and SO2 achieves consistent O3 values inside and outside the plumes. The achieved internal consistency indirectly validates the SO2 columns. In Fig. 16e, we show maps of DVCF-retrieved SO2 columns from a series of eight consecutive EPIC observations of the Raikoke plume in Fig. 16e, with the maximum SO2 value, the total SO2 mass, and the total area covered by elevated SO2 displayed in each snapshot. These results illustrate the high-cadence observing capability and high-quality SO2 measurements of EPIC.

6 Error analysis

We describe in this section how algorithm physics treatments and various sources contribute to the retrieval uncertainties and provide error estimates of the EPIC O3 and SO2 products.

6.1 General expression

The spectral measurements, represented by a column vector y of length m (the number of wavelength bands), are written explicitly with all the dependent parameters and possible errors and then expanded with respect to the linearization point (xL):

(31) y = ln I m = ln I TOA ( ω , ξ , b ) + ϵ m

(32) = ln I TOA ( ω L , ξ L , b L ) + k ω ( ω - ω L ) + k ξ ( ξ - ξ L ) + k b ( b - b L ) + ϵ f + ϵ m ,

where the column vectors of length nl (the number of atmospheric layers), ω, and ξ respectively represent the actual vertical profiles for O3 and SO2, while ωL is the climatological O3 profile equal to XmL)+γ1e1L) and ξL the prescribed SO2 profile specified by a GDF layer with an integrated vertical column equal to ΞL. The profile weighting functions, kω=-lnITOAω|ω=ωL and kξ=-lnITOAξ|ξ=ξL, are m×nl matrices, with each of the rows equal to the product of absorber cross-sections at one spectral band and the layer AMFs (i.e., the mean photon path lengths through the atmospheric layers). Likewise, b and bL are respectively a set of the exact forward model parameters and those used in the linearization, with the corresponding sensitivity matrix kb=-lnITOAb|b=bL. These forward model parameters may include the spectral-dependent MLER parameters, the ground surface and the OCP cloud pressures, the atmospheric temperature profile, the absorption cross-sections of O3 and SO2, and the parameters that specify the ISRFs of the spectral bands. The column vector ϵf is the forward modeling errors of the spectral bands, such as the approximate radiative transfer through Earth's spherical atmosphere and the incomplete accounting for RRS contributions. The last term ϵm is a column vector representing the spectral radiance errors of the instrument, including random noises and radiometric calibration biases.

Using the definition Δx=x-xL and putting Eq. (32) into residual ΔyL, Eq. (25) is rewritten as


where Aω=Gkω and Aξ=Gkξ are the averaging kernels (AKs) for O3 and SO2, respectively. Equation (34) describes how various error sources, from mismatches in absorber profiles to errors in model and measurement, propagate into the final result (x). The rows associated with the O3 and SO2 columns can be extracted from the vector equation (Eq. 34) and written as


after subtracting ΩT−ΩL and ΞT−ΞL from the row equations, respectively. Here ΩT and ΞT are the true O3 and SO2 columns, integrated from the corresponding true O3 (ω) and SO2 (ξ) profiles. AΩ and GΩ are the row vectors associated with the retrieved O3 column Ω from the corresponding matrices Aω and G. Analogously, AΞ and GΞ are the row vectors related to the retrieved SO2 column Ξ taken from the matrices Aξ and G, respectively. The constant row vector 1 contains the value 1 for all its elements. Thus, its dot product with a vertical profile (a column vector) is equivalent to the summation of all the individual layer amounts, yielding the total column. The column vector ϵΩ represents the total error combined from various sources impacting the total O3 accuracy, including errors in model parameters kb(bbL), forward modeling ϵf, spectral measurements ϵm, and the other absorber kξ(ξξL). Similarly, ϵΞ represents the combined total error affecting total SO2 accuracy with the other absorber term being replaced by kω(ωωL).

Retrieval errors can be characterized using Eqs. (35) and (36), provided errors from various sources are sufficiently small that forward modeling responds linearly to these deviations. However, substantial retrieval errors usually result from simplified physics treatments, which constrain the forward model to be radiative transfer in a molecular atmosphere over Lambertian surfaces. These errors may be called the AMF errors because the simplified physics treatments cannot, in general, reproduce the paths of photons through the observed atmosphere, even though they enable radiance matching between measurement and modeling. The deviations of mean paths lead to retrieval errors in O3 and SO2 because the interpretation of measured radiance through radiance matching requires accurate modeling of mean photon paths (i.e., the AMFs). The retrieval errors from the simplified physics treatment can be estimated using closed-loop tests (i.e., realistic forward modelings and then inverse retrieval with simplified physics treatments). Next, we provide uncertainty estimates of O3 and SO2 retrievals contributed from various error sources and simplified physics treatments.

6.2 Uncertainty estimates

6.2.1 Measurement errors

Errors in EPIC spectral measurements contribute to uncertainties in retrieved O3 and SO2 columns (Ω and Ξ). Taking the terms associated with radiance errors from Eqs. (35) and (36), retrieval errors are written as


These equations specify how measurement errors (ϵm) of the O3-sensitive bands and the MLER parameter errors (ΔR) due to the measurement errors in the weak absorption bands propagate into retrieved vertical columns.

Biases in radiance measurements lead to systematic errors in retrieved vertical columns. While the actual radiance biases are unknown, they are likely less than 1 % for EPIC UV bands. For radiance biases within ±1 %, the systematic O3 and SO2 column errors are within ±15 and ±8 DU, respectively, as estimated from Eqs. (37) and (38). These retrieval column errors are primarily controlled by the relative differences of spectral errors without significant dependence on the column amounts or surface reflectance. The retrieval biases vary with observing conditions given the same percentage radiance errors due to gain matrices (GΩ and GΞ) significantly depending on viewing and illumination angles.

In addition to systematic errors, radiance measurement noises add random errors to the retrieved columns. Retrieval errors due to random radiance noises (specified with normal distributions) are unbiased, with mean values, μ(ΔΩ) and μ(ΔΞ), close to zero and standard deviations, σ(ΔΩ) and σ(ΔΞ), proportional to standard deviations of radiance noises. The signal-to-noise ratios for EPIC UV bands are 290:1 (Herman et al.2018), equivalent to a noise level (standard deviation) of 0.345 % (1/290). This level is consistent with high-frequency radiance fluctuations (with standard deviations equal to 0.373 %, 0.354 %, 0.354 %, and 0.368 % for B1 to B4, respectively) within cloud-free scenes observed by EPIC. With a setting of equal standard deviations for the four UV bands (i.e., ϵm={0.345%,0.345%,0.345%,0.345%}), the estimated column O3 noise level is σ(ΔΩ)≃3.2 DU at low viewing zenith angles, decreases gradually with higher zenith angles, reaches a minimum of ∼1.5 DU at 75, and then rebounds quickly with further increases in zenith angles (see Fig. 17). The noise level of SO2 columns, σ(ΔΞ), exhibits a similar angular dependence as shown in Fig. 17, primarily following the angular variation of the gain matrix GΞ. These angular-dependent column noises are insensitive to the column amounts or the surface reflectance.

Figure 17Noise levels, i.e., standard deviations (σ) of O3 and SO2 errors (Δ), contributed from the random noises in EPIC spectral measurements. The SO2 noise estimate is for a layer at an altitude 11 km above sea level.


6.2.2 Model parameter errors

Retrieval errors due to uncertainties in model parameters, including molecular absorption cross-sections (σ) and atmospheric temperature profiles (T), are estimated as


where (mzΔσO3,SO2) and (mzσO3,SO2T) are m×nl matrices with elements {mzj(λi)ΔσO3,SO2(λi,Tj),i=1m,j=1nl} and {mzj(λi)σO3,SO2(λi,T)T|Tj,i=1m,j=1nl}, respectively. Here mz (see Eq. 7) is the mean photon path length through a layer at altitude z, ΔσO3,SO2 represents errors in O3 or SO2 cross-sections, and ΔT represents errors in atmospheric temperature profiles.

The BDM O3 cross-sections (Daumont et al.1992; Brion et al.1993; Malicet et al.1995, commonly abbreviated as BDM) and the BW SO2 cross-sections (Birk and Wagner2018) are used in O3 and SO2 retrievals from EPIC. These baseline cross-sections contain errors, which are not known quantitatively but can be estimated by comparing with alternative cross-sections. Specifically, the BW O3 cross-sections (Birk and Wagner2021) and the SO2 absorption cross-sections of Bogumil et al. (2003) are the alternatives that can replace the baselines for EPIC retrievals. The cross-section errors, ΔσO3 in Eq. (39) and ΔσSO2 in Eq. (40), are estimated based on the differences between the alternatives and baselines, showing that alternative O3 and SO2 cross-sections are slightly (about 0.1 % to 1.1 %) lower than the corresponding baselines at EPIC B1 and B2. These biases in cross-sections result in O3 column biases between 0.5 % and 2 % as well as SO2 column biases between 1 % and 2 % depending on the effective cross-section differences. The temperature dependence of the BW O3 cross-sections behaves quite differently from the BDM (Bak et al.2020), especially at EPIC B2, contributing to the high ends of O3 biases (>1.0 %), which occur predominantly at high (viewing or solar) zenith angles when O3 retrieval becomes more sensitive to EPIC B2 radiance.

Both O3 and SO2 cross-sections are significantly dependent on temperatures. Thus, accurate temperature profiles are needed to determine atmospheric absorption properties for modeling of measured radiances. As mentioned in Sect. 3, MERRA-2 climatological temperature profiles (Yang and Liu2019) are used for retrievals from EPIC. Actual temperature profiles differ from the climatological profiles. Over a short period (e.g., a day), the spatial distribution of these differences is not random, leading to retrieval errors that are unevenly distributed spatially. However, actual temperature profiles are normally distributed around the climatological mean over a long period (e.g., a month) for a location. Therefore, temperature profile mismatches add random components, which average to zero over a long time, to the total errors. The variances of these random errors are proportional to the layer-column-weighted temperature error variances. Estimated from the variances of temperature profiles (Yang and Liu2019), the random components, σ(ΔΩ), are ∼0.3 % in the tropics, increase to ∼0.7 % in the midlatitudes, and reach ∼1 % at high latitudes. Similarly, random errors, σ(ΔΞ), are ∼0.8 % in the tropics, ∼1.7 % in the midlatitudes, and ∼3.5 % at high latitudes.

6.2.3 Forward modeling errors

The MLER treatment adopted for the retrieval algorithm allows the use of the vector radiative transfer code, TOMRAD, as the forward model to simulate measured radiances and weighting functions. TOMRAD implements Dave's iterative solution (Dave1964) with pseudo-spherical approximation (Caudill et al.1997) to the problem of the transfer of solar radiation through a molecular atmosphere over a Lambertian surface. The forward modeling with TOMRAD is accurate for EPIC observations around the center of its hemispheric view, with radiance errors (ϵf) of all EPIC UV bands less than ±0.2 % for VZA <50. Note that for EPIC observations, each of its IFOVs has similar VZA and SZA (see Fig. 1) with differences VZA-SZA<±9. As EPIC observations move towards the edge, the pseudo-spherical model atmosphere deviates more from Earth's spherical atmosphere in accounting for atmospheric attenuation and multiple scattering, resulting in more significant errors in modeled radiances, whose maximum errors increase to about ±1 % at 75 VZA and about ±2 % at 85 VZA (Caudill et al.1997). RRS corrections are included in the forward modeling, and they are well within ±1 % for EPIC UV bands (see, e.g., Fig. 10). Incomplete RRS corrections are expected to add less than ±0.1 % to the forward modeling errors.

Unlike the calibration biases being insensitive to observing conditions and having no correlation among different bands, the forward modeling errors vary with absorber amounts and surface reflection, and they overestimate or underestimate similarly for all the UV bands depending on the viewing and illumination geometry. How these radiance errors propagate into the retrieved columns can be estimated using Eqs. (37) and (38), with error source terms replaced by ϵf and ΔR due to modeling errors in the long wavelength bands. With the radiance errors estimated above, these equations yield retrieval errors up to ±0.6 and ±0.3 DU when VZA <50, increasing to ±1.5 and ±1 DU at VZA =75, as well as ±5 and ±15 DU at VZA =85, respectively, for O3 and SO2 vertical column errors. These are systematic errors and vary between high and low biases spatially depending on observing conditions, especially the viewing and illumination geometry.

6.2.4 Profile errors

As described in Sect. 3, a column-dependent O3 profile, whose tropospheric integration matches the climatological tropospheric column, is used to specify the vertical distribution of a retrieved total O3 column. This retrieved profile (ωL), which represents the likely vertical distribution of the retrieved O3 vertical column, invariably differs from the actual profile (ω). The O3 error (ΔΩ) due to a profile errors (ωωL) can be quantified using the first term on the rhs of Eq. (35), which is regulated by the retrieval AK (AΩ). Examples of AKs for EPIC total O3 retrievals are shown in Fig. 18a and c, illustrating how AΩ changes with viewing geometry. For low VZAs (<55), O3 AKs are close to 1 above the upper troposphere, and therefore profile mismatches in this altitude region result in insignificant retrieval errors. However, profile mismatches produce sizable retrieval errors for high VZAs. In the troposphere, O3 AKs change with surface reflectance in addition to angular dependence. Under cloud-free conditions, O3 AKs drop quickly towards low-reflectivity surfaces, more so at high zenith angles (see, e.g., Fig. 18a). Above highly reflective surfaces (e.g., snow, ice, or bright clouds), O3 AKs increase drastically (see Fig. 18c), indicating enhanced sensitivity to tropospheric profile, especially at low zenith angles. Evidently, O3 errors due to profile mismatches primarily come from the troposphere for low zenith angles for both low- and high-reflectivity surfaces, though stratospheric contributions increase substantially with higher zenith angles, more significantly for a low-reflectivity surface. In general, errors due to the profile shape are reduced for high-reflectivity surfaces.

Figure 18Examples of EPIC total O3 AKs (a, c) and SO2 AKs (b, d) as functions of geometric altitude (z) above seal level for several VZAs (θv). These AKs are calculated for a molecular atmosphere at midlatitude with 275 DU total O3 and 30 DU of SO2 in a layer (1.5 km thick) at 11 km altitude over low-reflectance (a, b) and high-reflectance (c, d) Lambertian surfaces. Observing conditions are listed at the top of the panels.


Over a short period (e.g., 1 d), O3 errors due to profile mismatches are local biases (reductions or enhancements) that vary smoothly with location. However, they are random errors since mismatches (ωωL) are normally distributed around their near-zero means over a long period (e.g., 1 month). The variances of O3 errors (ΔΩ) can be written as

(41) Var Δ Ω = E A Ω - 1 ω - ω L 2 = A Ω - 1 E ω - ω L ω - ω L T A Ω - 1 T = A Ω - 1 S n l A Ω - 1 T ,

where the expected values (i.e., the statistical means), E(ω-ωL)(ω-ωL)T, are O3 profile covariance matrices Snl, which depend on total columns (Ω), season, and latitude. This random component is estimated as a function of VZA using the column-dependent Snl from the M2TCO3 climatology (Yang and Liu2019). Figure 19 shows that the standard deviation of this error component increases gradually with higher VZA, from 1 % at nadir to 1.7 % at 75, then rapidly with further elevated VZA.

Figure 19The standard deviation (σ) of O3 errors due to profile mismatches as a function of the viewing zenith angle (θv), estimated using the M2TCO3 profile covariance matrix for the December midlatitude zone (40–50 N).


Retrieval of SO2 requires knowledge of the altitude at the center of the volcanic plume, which can be represented by a narrow (e.g., a width of 1.5 km) GDF. Error in the plume altitude leads to SO2 retrieval error, which can be estimated using the retrieval AK. Figure 18b and d show sample AKs and their variations with VZA for an SO2 plume center at 11 km altitude. The values of these AKs are equal to 1 at 11 km, meaning no retrieval error when the altitude used in the retrieval is equal to the actual plume altitude. Here, we examine the case of a low-reflectivity surface. The retrieved SO2 column overestimates (underestimates) the actual column when the plume altitude is higher (lower) than the assumed altitude (11 km). At low VZAs (<35), AK values are close to 1 above the assumed altitude (11 km), indicating small (less than a few percents) errors for plumes at higher altitudes. Overestimation (up to 150 %) increases quickly with larger VZAs when the plume is at higher altitudes. Underestimation is more severe with a lower altitude of the plume: 10 % to 20 % per 1 km lower than the assumed altitude.

6.2.5 Errors from Lambertian treatment of natural surfaces

Reflections from surfaces are anisotropic but treated as isotropic. To estimate errors in O3 and SO2 columns due to this simplification, we performed DVCF retrieval from simulated radiances. First, TOA radiances of the four EPIC UV bands are modeled using a state-of-art radiative transfer model, VLIDORT (Spurr2006), for a molecular atmosphere with various O3 and SO2 profiles over a surface characterized by an anisotropic BRDF. Next, GLERs are determined at the long wavelength bands (B3 and B4) and then linearly extrapolated to the short wavelength bands (B1 and B2). Finally, retrieved O3 and SO2 columns from simulated B1 and B2 radiances using the extrapolated GLERs are compared with the column settings of the forward modeling to quantify retrieval errors. Examples of O3 and SO2 errors determined this way are shown in Fig. 20 for observing conditions described in the Fig. 4 caption. In the closed-loop testing, surface reflection is specified by the Cox–Munk BRDF, which is highly anisotropic, more so than the land surface BRDFs that are well-characterized by the combinations of Lambertian, Ross, and Li kernels (Lucht et al.2000; Schaaf et al.2011). Hence, the Cox–Munk BRDF selection provides error ranges due to the Lambertian treatment of surface reflections. Closed-loop tests are performed for a wide range of viewing illumination geometries and vertical distributions of O3 and SO2. Test results (see, e.g., Fig. 20) show that errors in total O3 are mostly within ±1 DU, while SO2 errors are within ±5 % for SO2 layers above 5 km, and decrease (increase) with higher (lower) layer altitudes. As shown in Fig. 5, the AMF errors due to Lambertian treatment occur below 20 km altitude. Consequently, a small fraction of the O3 profile is affected by this approximation. Thus, O3 errors are proportional to the tropospheric columns but are insensitive to the total column amounts. Since a vast majority of volcanic SO2 clouds are located below 20 km altitude, SO2 errors are proportional to the total SO2 columns. Higher SO2 clouds are not affected by this treatment.

Figure 20Errors in retrieved O3 and SO2 due to Lambertian surface treatment of an anisotropic surface. (a) O3 errors in DU for a midlatitude O3 profile with total columns of 275, 375, and 475 DU. (b) SO2 errors in percent for an SO2 layer at altitudes of 5, 7, and 11 km above sea level. See the Fig. 4 caption for the specification of surface BRDF and viewing and illumination geometry.


6.2.6 Errors from MLER treatment of clouds and aerosols

The MLER model is adopted to treat atmospheric particles, including clouds and aerosols, which predominantly reside in the lower troposphere. The modeled light paths (especially in the troposphere) based on this treatment differ significantly from those for light transfer through the particle-laden atmosphere (see, e.g., Fig. 9). The retrieval errors due to this simplification are again estimated using closed-loop testing. First, TOA radiances of the EPIC UV bands are simulated using VLIDORT for particle-laden atmospheres with various O3 and SO2 profiles over Lambertian surfaces of different reflectivities. Then inversion from the simulated radiances with the MLER treatment permits the identification of conditions under which retrieval errors are significant.


The error in total O3 due to the MLER treatment of a low-lying (below 10 km) cloud is mostly within ±2 DU (e.g., Fig. 21a). This O3 error decreases slightly with a lower cloud altitude (or higher cloud pressure) but is insensitive to the cloud fraction or the total O3 column. In other words, the MLER treatment does not contribute to large uncertainty in the retrieved total O3 column, provided that an accurate OCP for the cloud is used for the MLER cloud surface. However, OCP has some uncertainty, contributing to additional uncertainty in the O3 column: a low (high) bias in OCP results in a positive (negative) error in total O3; quantitatively, ±100 hPa causes about ∓4 DU (see Fig. 21a). The OCP uncertainty is estimated to be within ±50 hPa, thus contributing ∓2 DU to the total O3 uncertainty. Combining O3 uncertainties due to the OCP error and the MLER treatment yields ±4 DU uncertainty in total O3 under cloudy conditions.

Figure 21Errors in retrieved O3 and SO2 due to MLER treatment of clouds, which are represented by 1.5 km thick C1 particle layers with an optical thickness τ=15 at 340 nm. (a) O3 errors in DU for correct (pc=545 hPa) and biased (pc=545±100 hPa) cloud OCPs. (b) SO2 errors in percent for SO2 layers at three altitudes (5, 7, and 11 km) above a layer of cloud (at 3 km altitude). See the Fig. 4 caption for the specification of viewing and illumination geometry.


The error in the total SO2 column due to the MLER cloud treatment is within ±2 % when the SO2 layer in the troposphere is well above the underlying cloud. This SO2 error increases with a smaller separation between the SO2 layer and the cloud, reaching ±15 % when the SO2 layer is just above the cloud. These characteristics of SO2 error are illustrated in Fig. 21b. In contrast to the MLER treatment O3 error, which is insensitive to the total column, this SO2 error is proportional to the total SO2 column. When the SO2 layer is below or within the cloud, the uncertainty of SO2 quantification increases drastically. Depending on the relative distributions of SO2 and the cloud particles, the retrieved SO2 based on the MLER treatment can be a fraction of or a few times the actual column.


Besides clouds, the MLER treatment is applied to IFOVs contaminated by aerosols, which primarily reside in the troposphere and cover a significant portion of Earth's surface. These aerosols are suspended tiny (micron-scale) particles that scatter and possibly absorb sunlight. The frequently observed non-absorbing (or weakly absorbing) aerosols are sea salt and sulfate (SLF), and UV-absorbing aerosols are smoke (i.e., carbonaceous aerosols from biomass combustion, BIO), mineral dust (DST), and volcanic ash. Moderate and high positive AI values indicate the presence of UV-absorbing aerosols in an IFOV, while negative and slightly positive AI values indicate the presence of non-absorbing or weakly absorbing aerosols.

Closed-loop testing shows that MLER treatment of non-absorbing and weakly absorbing aerosols in the lower troposphere (<7 km) results in small (<±2 DU) errors in total O3 retrievals, provided that the proper OCP for the elevated cloud surface is used. This error range is nearly independent of the total O3 column or the aerosol loading.

The MLER treatment errors for UV-absorbing aerosols close to the surface (<1 km altitude) are mostly within ±1 DU, similar to the error range associated with the LER treatment of BRDF surfaces. For elevated UV-absorbing aerosols, the MLER treatment and the linear rgfc extrapolation scheme (see Sect. 2.3) result in a positive bias in the retrieved total O3 columns (see, e.g., Fig. 22). This O3 bias depends on the viewing illumination geometry and generally increases with stronger aerosol absorption (i.e., lower single-scattering albedo, ω), larger aerosol optical thickness, and higher altitude of aerosol layers. Regression analysis of results from closed-testing with many combinations of viewing illumination geometries, particle-laden atmospheres (with various optical properties, optical thicknesses, and vertical distributions), surface reflectivities, and O3 profiles reveals a positive correlation between column O3 error and the UV AI. Quantitatively, this relationship can be written as ΔO3=(1.5±1)×AI DU for AI values greater than 0.5 and less than 8. This relationship provides a rough estimate of O3 bias based on the observed AI. Typically, AI values fall between 1 and 4 with a median value of 1.5 for EPIC observations of UV-absorbing aerosols, corresponding to a mean O3 bias of about 3 DU for IFOVs contaminated with UV-absorbing aerosols. The MLER treatment sometimes fails when aerosol absorption is strong such that the derived LER becomes negative. In this case, the explicit aerosol treatment may be needed to reduce the retrieval uncertainty.

Figure 22Errors in retrieved O3 due to MLER treatment of two common UV-absorbing aerosols, (a) BIO (ω=0.921) and (b) DST (ω=0.900), with various optical thicknesses (τ = 0.25, 0.5, 1.0, and 2.0 at 340 nm) located at 5 km altitude. See the Fig. 4 caption for the specification of viewing and illumination geometry. The AIs associated with each observation scenario are shown in (c) and (d).


Figure 23Intercomparison of total O3 from EPIC and the ground-based Brewer spectrophotometers from July 2015 to April 2021 at 10 selected ground stations with high-cadence measurements: Alert (82.50 N), Eureka (79.99 N), Resolute (74.72 N), Churchill (58.75 N), Edmonton (53.55 N), Goose Bay (53.31 N), De Bilt (52.10 N), Thessaloniki (40.63 N), Paramaribo (5.806 N), and the South Pole (−89.99 N). EPIC and Brewer coincident pairs are used in the plots, and data with VZA 70 N only are included in the difference statistics.


EPIC's high-cadence observation has more chances to view volcanic clouds during or soon after eruptions. These young volcanic clouds contain mixtures of ash particles and water or ice clouds, as eruptions inject ash and gases (including SO2) into the atmosphere. Since ash particles strongly absorb UV, the MLER treatment of volcanic plumes leads to huge uncertainties in the retrieved SO2 columns, which are often greatly overestimated or underestimated depending on the relative distributions between SO2 and ash particles. An explicit treatment of volcanic ash is needed for accurate retrieval of SO2 when ash particles are co-located with or slightly separated from the gas.

Figure 24Comparison of synoptic EPIC O3 with MERRA-2 assimilated O3: time series of mean daily differences and standard deviations for EPIC observations with VZA ≤70.


6.3 Error summary

Uncertainty estimates (Sect. 6.2) have detailed various contributions to systematic and random errors in retrieved O3 and SO2 vertical columns. In this summary, these error types are separately combined to estimate their total systematic and random errors, respectively.

First, random errors from various sources, including measurement noise and errors in temperature profiles and O3 vertical profiles, are combined to estimate total random errors. The random O3 error (characterized by its standard deviation) is 1.5 % at low VZAs (<45), increasing to 2.0 % at 75 for IFOVs without clouds and aerosols. Similarly for IFOVs with clouds or non-absorbing aerosols, the random error in the total O3 column is 1.8 % for low VZAs (<45) and 2.5 % at 75. Random SO2 error has two terms: one is independent of the SO2 column (see Sect. 6.2.1), but the other is proportional to this column (see Sect. 6.2.2). Hence, they are not combined.

Figure 25EPIC and OMPS observations of volcanic SO2 plumes on 17 June 2018 from the eruption of Fernandina volcano (Δ) in the Galapagos Islands. This eruption injected a significant amount of SO2 into the troposphere at about 3.5 km above sea level. The mass loading of an SO2 plume is obtained by summing the SO2 masses of all IFOVs with SO2 vertical columns ≥1 DU. The lower right panel plots the EPIC and OMPS SO2 masses vs. the observation time (UTC).

Next, the radiometric biases, forward radiance modeling errors, and errors in molecular cross-sections are the primary contributions to the systematic errors in the retrieved vertical columns. The possible ranges due to these sources are provided earlier in this section, but their actual contributions are unknown. We obtain the total systematic error by combining the likely ranges (i.e., ∼half of the possible ranges) of these contributions and estimate the bias in the total O3 column to be ±2 % for VZA 70 and ±3.5 % for VZA 85; the bias in the total SO2 column high (>10 km) in the atmosphere is estimated to be ±1.5 % for VZA 70 and ±5 % for VZA 85. Some biases depend on geophysical conditions. For instance, O3 has a positive bias for IFOVs with UV-absorbing aerosols, increasing with a higher AI (see Sect. 6.2.6). Underestimation (overestimation) of an SO2 column occurs when its layer altitude used in the retrieval is higher (lower) than the actual altitude (see Sect. 6.2.3).

7 Validation of EPIC O3 and comparison of SO2

7.1 O3 validation

We validate the DVCF O3 retrievals from EPIC using ground-based Brewer spectrophotometer measurements and the assimilated O3 product from MERRA-2, the Modern-Era Retrospective Analysis for Research and Applications Version 2 (Gelaro et al.2017).

We compare EPIC total O3 columns with the Brewer O3 data at 10 selected ground stations with high-cadence measurements distributed in five latitude zones. At each of these selected stations, a Brewer spectrometer makes a measurement every few (∼10) minutes during the daylight hours each day, thus providing total vertical O3 columns that are coincident (within ±15 min) with EPIC observations at these stations. For intercomparison, Brewer O3 data are interpolated to the times when EPIC observes these locations. Coincident O3 data from these two independent sources are displayed in the upper panels of Fig. 23a–e, their differences in the lower panels of Fig. 23a–e, and the EPIC vs. Brewer scatter plots in the right panels of Fig. 23a–e. We include coincident data with VZA ≤70 only for statistical analysis due to EPIC–Brewer IFOV differences that usually increase with slant path lengths and EPIC's footprint sizes as well as due to calibration biases at large VZAs or SZAs. The mean difference and standard deviations in percent are displayed in the difference plots, while those in DU and the correlation coefficients are in the scatter plots. Time series of O3 difference (see lower panels of Fig. 23a–e between EPIC and Brewer) are highly stable with similar moving averages and standard deviations from June 2015 to April 2021, showing that EPIC O3 is consistent over time, without noticeable drift. The correlations between EPIC and Brewer are very high with correlation coefficients R≥0.96 for most stations (except for the Paramaribo station near the Equator, where R=0.87), demonstrating that EPIC captures O3 variability accurately. EPIC O3 agrees with the Brewer measurements to better than 1 % with standard deviations of differences less than 3.5 % for all the ground stations, validating the high accuracy of EPIC total O3.

From October 2004, MERRA-2 O3 field is assimilated from Aura MLS and OMI and provides highly realistic global distributions of O3 in the stratosphere and upper troposphere while inheriting the uncertainty characteristics of its sources (Stajner et al.2008; Wargan et al.2015; Davis et al.2017). We compare the MERRA-2 synoptic O3 field with the EPIC hemispheric view for the same observation time to access EPIC's capability of capturing the realistic O3 distribution. For instance, strikingly similar O3 spatial distributions are observed in EPIC measurements (Fig. 15c) and the MERRA-2 assimilation (Fig. 15d), with agreement at -0.20±2.52 % (or -0.35±5.6 DU; Fig. 15h). We expand this synoptic comparison to each EPIC hemispheric view obtained from July 2015 to August 2021 and plot in Fig. 24 the time series of daily statistics. This time series shows that nearly the same level of agreement is achieved for the entire period, with a mean bias and standard deviation of 0.44±2.19 % (or 1.16±6.34 DU). Considering the mean bias (about −1.2 %) of MERRA-2 total O3 (Wargan et al.2017), we estimate the accuracy of EPIC total O3 to be -0.76±2.19 %.

7.2 SO2 comparison

Volcanic eruptions occur sporadically and without warning, but EPIC on DSCOVR, from the unique L1 vantage point, usually provides multiple daily observations of volcanic SO2 and ash clouds once injected into the atmosphere. In contrast, ground-based instruments rarely detect volcanic clouds unless they drift over one in operation. We thus rely on polar-orbiting instruments, which may observe a volcanic cloud once (or more at high latitude) per day to provide validation measurements.

The OMPS-NM on SNPP provides high-quality hyperspectral measurements in the UV, from which highly accurate retrievals of O3 and SO2 are achieved using the DVCF algorithm. The DVCF algorithm can apply to both discrete spectral measurements (e.g., TOMS and EPIC) and hyperspectral ones (e.g., OMI and OMPS-NM). The main difference is that more information can be extracted from hyperspectral measurements to improve the accuracy and precision of the retrieved geophysical parameters. For instance, the altitude of an SO2 layer can be determined in addition to its amount simultaneously using the DVCF algorithm (Yang et al.2010). Having the altitude information significantly improves the accuracy of SO2 quantification because the SO2 measurement sensitivity varies strongly with its altitude. Thus, DVCF height-resolved SO2 retrievals from hyperspectral instruments, such as OMI and OMPS-NM, provide the most accurate quantification of SO2 vertical columns. To validate DVCF SO2 retrievals from EPIC, we compare the SO2 mass loading of a volcanic plume integrated from EPIC observations with SNPP OMPS-NM for the same event.

Figure 25 compares the DVCF retrievals of the volcanic plume from the explosive eruption of Fernandina volcano in the Galapagos Islands on 17 June 2018. Seven plume exposures about 65 min apart are taken by EPIC on this day. Soon after the fifth EPIC exposure, the OMPS observed this plume for the first time. For the exposures at 10 min apart, both instruments estimate the mass loading at 71 kt, validating the EPIC SO2 result.

The lower right panel of Fig. 25 plots the plume mass vs. the observation time, showing the mass loading peaks near local noontime. The observed mass change results from the continuing emission from the volcano, the conversion of SO2 into sulfate, and the changing measurement sensitivity with viewing illumination conditions since low SO2 columns may be missed at large (VZA, SZA, or both) angles due to lower sensitivity. EPIC's high-cadence observations allow better identification of the peak loading of volcanic SO2 plume, thus usually providing more accurate estimates of the lower bound of SO2 emission compared to polar-orbiting instruments.

We have conducted many mass loading comparisons between EPIC and OMPS and found that the agreements are usually within 20 %. These findings indicate that DVCF SO2 retrieval from EPIC provides better than 20 % (an estimate of the upper error bound) accuracy in total mass for eruptions with greater than 50 kt emissions.

8 Conclusions

We present the algorithm for making the EPIC O3SO2AI product in this paper. This algorithm is based on the DVCF algorithm developed for retrieving trace gases, including O3, SO2, and NO2, from Aura OMI and Suomi NPP as well as NOAA-20 OMPS. Algorithm advances, including the improved O3 profile representation and the regulated direct fitting inversion technique, improve the accuracy of O3 and SO2 from the multi-channel measurements of DSCOVR EPIC. The theoretical basis of the SOE approach, introduced to reduce retrieval artifacts due to EPIC's band-to-band misregistration, can be exploited for other applications, such as the separation of a spatially smooth data field (e.g., stratospheric O3) from that (e.g., tropospheric O3) with higher spatial variations.

A thorough error analysis is provided to quantify O3 and SO2 retrieval uncertainties due to various error sources and simplified algorithm physics treatments. Error analysis findings indicate that the MLER treatment of UV-absorbing aerosols leads to significant uncertainties in retrieved O3 and SO2 columns. Future improvements may include explicit aerosol treatment or other schemes for radiance or product corrections. The GLER treatment of anisotropic surface reflections introduces small errors in the retrieved total O3 and SO2 columns, primarily because surface reflection is a minor component of measured radiance in the UV. However, this GLER treatment does not generally provide a more accurate tropospheric AMF. Hence, explicit BRDF treatment of surface reflection is needed for accurate retrievals of tropospheric gases.

The EPIC total O3 columns are validated against coincident ground-based Brewer measurements and compared with coincident O3 data from MERRA-2 assimilation. The findings show that EPIC total O3 is highly accurate, realistically capturing the short-term O3 variability while maintaining long-term consistency over the entire record. The EPIC SO2 loadings of volcanic plumes are evaluated against those from hyperspectral measurements of the same eruptions, showing that EPIC provides accurate SO2 quantifications from large volcanic eruptions. EPIC's high-cadence observations allow better identification of the peak loading of volcanic SO2 plumes compared to polar-orbiting instruments.

Data availability

The EPIC product, O3SO2AI (NASA/LARC/SD/ASDC2018), contains scene reflectivity, aerosol index (AI), total vertical columns of ozone (O3), and vertical sulfur dioxide (SO2) columns when volcanic clouds are detected in the EPIC field of view. This product is available from the NASA Langley Atmospheric Science Data Center (ASDC), accessible at this link: (last access: 1 October 2022).

The BW O3 and SO2 cross-sections (Birk and Wagner2021, 2018) are available at the Zenodo repository with DOIs and, respectively.

Author contributions

XH (the lead author) contributed to the development of computer codes, conducted data analysis and product validations, and prepared the initial draft of the paper. KY conceived the research, designed the algorithm, contributed to algorithm implementation, data analysis, and product validations, and contributed to writing the article.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


We are grateful to the DSCOVR science team for providing EPIC L1 data and to the NASA Center for Climate Simulations for providing resources for data processing. We would like to thank the editor M. Van Roozendael and two anonymous reviewers for their helpful comments and suggestions that improved the paper.

Financial support

This research has been supported by the NASA DSCOVR Program managed by Richard S. Eckman (grant nos. 80NSSC22K0502, 80NSSC19K0770, and NNX15AB10G).

Review statement

This paper was edited by Michel Van Roozendael and reviewed by two anonymous referees.


Ahmad, Z., Bhartia, P. K., and Krotkov, N. A.: Spectral properties of backscattered UV radiation in cloudy atmospheres, J. Geophys. Res.-Atmos., 109, D01201,, 2004. a, b

Bak, J., Liu, X., Wei, J. C., Pan, L. L., Chance, K., and Kim, J. H.: Improvement of OMI ozone profile retrievals in the upper troposphere and lower stratosphere by the use of a tropopause-based ozone profile climatology, Atmos. Meas. Tech., 6, 2239–2254,, 2013. a

Bak, J., Liu, X., Birk, M., Wagner, G., Gordon, I. E., and Chance, K.: Impact of using a new ultraviolet ozone absorption cross-section dataset on OMI ozone profile retrievals, Atmos. Meas. Tech., 13, 5845–5854,, 2020. a

Bhartia, P. K. and Wellemeyer, C. G.: TOMS-V8 Total O3 Algorithm, in: OMI Algorithm Theoretical Basis Document, 2nd edn., vol. II, edited by: Bhartia, P. K., NASA Goddard Space Flight Center, Greenbelt, Maryland, USA, 15–32, (last access: 1 October 2022), 2002. a, b, c, d, e

Birk, M. and Wagner, G.: ESA SEOM-IAS – Measurement and ACS database SO2 UV region (Version 1), Zenodo [data set],, 2018. a, b

Birk, M. and Wagner, G.: ESA SEOM-IAS – Measurement and ACS database O3 UV region (Version 2), Zenodo [data set],, 2021. a, b

Blank, K.: EPIC Geolocation and Color Imagery Algorithm Revision 6, NASA Goddard Space Flight Center, Greenbelt, Maryland, USA,, 2019. a

Blank, K., Huang, L.-K., Herman, J., and Marshak, A.: Earth Polychromatic Imaging Camera Geolocation; Strategies to Reduce Uncertainty, Front. Remote Sens., 2, ISSN: 2673-6187,, 2021. 

Bogumil, K., Orphal, J., Homann, T., Voigt, S., Spietz, P., Fleischmann, O., Vogel, A., Hartmann, M., Kromminga, H., Bovensmann, H., Frerick, J., and Burrows, J.: Measurements of molecular absorption spectra with the SCIAMACHY pre-flight model: instrument characterization and reference data for atmospheric remote-sensing in the 230–2380 nm region, J. Photoch. Photobio. A, 157, 167–184,, 2003. a

Bosilovich, M., Akella, S., Coy, L., Cullather, R., Draper, C., Gelaro, R., Kovach, R., Liu, Q., Molod, A., Norris, P., Wargan, K., Chao, W., Reichle, R., Takacs, L., Vikhliaev, Y., Bloom, S., Collow, A., Firth, S., Labow, G., Partyka, G., Pawson, S., Reale, O., Schubert, S. D., and Suarez, M.: MERRA-2: Initial Evaluation of the Climate, in: NASA Technical Report Series on Global Modeling and Data Assimilation NASA/TM–2015-104606, vol. 43, edited by: Koster, R. D., Goddard Space Flight Center, Greenbelt, Maryland, USA, p. 139, (last access: 1 October 2022), 2015. a

Bovensmann, H., Burrows, J. P., Buchwitz, M., Frerick, J., Noël, S., Rozanov, V. V., Chance, K. V., and Goede, A. P. H.: SCIAMACHY: Mission Objectives and Measurement Modes, J. Atmos. Sci., 56, 127–150,<0127:SMOAMM>2.0.CO;2, 1999. a

Brennan, B. and Bandeen, W. R.: Anisotropic Reflectance Characteristics of Natural Earth Surfaces, Appl. Optics, 9, 405,, 1970. a

Brion, J., Chakir, A., Daumont, D., Malicet, J., and Parisse, C.: High-resolution laboratory absorption cross section of O3. Temperature effect, Chem. Phys. Lett., 213, 610–612,, 1993. a

Brodzik, M. J. and Stewart, J. S.: Near-Real-Time SSM/I-SSMIS EASE-Grid Daily Global Ice Concentration and Snow Extent, Version 3, NASA National Snow and Ice Data Center Distributed Active Archive Center, Boulder, Colorado, USA,, 2021. 

Callies, J., Corpaccioli, E., Eisinger, M., Hahne, A., and Lefebvre, A.: GOME-2 – Metop's Second-Generation Sensor for Operational Ozone Monitoring, ESA Bulletin, (last acccess: 1 October 2022), 2000. a

Caudill, T. R., Flittner, D. E., Herman, B. M., Torres, O., and McPeters, R. D.: Evaluation of the pseudo-spherical approximation for backscattered ultraviolet radiances and ozone retrieval, J. Geophys. Res.-Atmos., 102, 3881–3890,, 1997. a, b

Cede, A., Kang Huang, L., McCauley, G., Herman, J., Blank, K., Kowalewski, M., and Marshak, A.: Raw EPIC Data Calibration, Front. Remote Sens., 2, 1–18,, 2021. 

Chance, K. V.: Spectroscopic Measurements of Tropospheric Composition from Satellite Measurements in the Ultraviolet and Visible: Steps Toward Continuous Pollution Monitoring from Space, in: Remote Sensing of the Atmosphere for Environmental Security, edited by: Perrin, A., Ben Sari-Zizi, N., and Demaison, J., Springer Netherlands, Dordrecht, 1–25,, 2006. a

Chance, K. V. and Spurr, R. J. D.: Ring effect studies: Rayleigh scattering, including molecular parameters for rotational Raman scattering, and the Fraunhofer spectrum, Appl. Optics, 36, 5224,, 1997. a

Coulson, K. L. and Reynolds, D. W.: The Spectral Reflectance of Natural Surfaces, J. Appl. Meteorol., 10, 1285–1295, (last access: 1 October 2022), 1971. a

Cox, C. and Munk, W.: Measurement of the Roughness of the Sea Surface from Photographs of the Sun's Glitter, J. Opt. Soc. Am., 44, 838–850,, 1954a. a

Cox, C. and Munk, W.: Statistics of the sea surface derived from Sun glitter, J. Mar. Res., 13, 198–227, 1954b. a

Daumont, D., Brion, J., Charbonnier, J., and Malicet, J.: Ozone UV spectroscopy I: Absorption cross-sections at room temperature, J. Atmos. Chem., 15, 145–155,, 1992. a

Dave, J. V.: Meaning of Successive Iteration of the Auxiliary Equation in the Theory of Radiative Transfer, Astrophys. J., 140, 1292,, 1964. a, b, c, d, e

Davis, S. M., Hegglin, M. I., Fujiwara, M., Dragani, R., Harada, Y., Kobayashi, C., Long, C., Manney, G. L., Nash, E. R., Potter, G. L., Tegtmeier, S., Wang, T., Wargan, K., and Wright, J. S.: Assessment of upper tropospheric and stratospheric water vapor and ozone in reanalyses as part of S-RIP, Atmos. Chem. Phys., 17, 12743–12778,, 2017. a

Deirmendjian, D.: Electromagnetic scattering on spherical polydispersions, American Elsevier Pub. Co., New York, ISBN: 0444000380, (last access: 1 October 2022), 1969. a

Doda, D. D. and Green, A. E. S.: Surface reflectance measurements in the UV from an airborne platform Part 1, Appl. Optics, 19, 2140,, 1980. a

Doda, D. D. and Green, A. E. S.: Surface reflectance measurements in the ultraviolet from an airborne platform – Part 2, Appl. Optics, 20, 636,, 1981. a

Eck, T. F., Bhartia, P. K., Hwang, P. H., and Stowe, L. L.: Reflectivity of Earth's surface and clouds in ultraviolet from satellite observations, J. Geophys. Res.-Atmos., 92, 4287,, 1987. a

Edlén, B.: The Refractive Index of Air, Metrologia, 2, 71–80,, 1966. a

Eskes, H. J., van der A, R. J., Brinksma, E. J., Veefkind, J. P., de Haan, J. F., and Valks, P. J. M.: Retrieval and validation of ozone columns derived from measurements of SCIAMACHY on Envisat, Atmos. Chem. Phys. Discuss., 5, 4429–4475,, 2005. a

Feister, U. and Grewe, R.: Spectral albedo measurements in the UV and visible region over different types of surfaces, Photochem. Photobiol., 62, 736–744,, 1995. a

Flynn, L. E., Long, C., Wu, X., Evans, R., Beck, C. T., Petropavlovskikh, I., McConville, G., Yu, W., Zhang, Z., Niu, J., Beach, E., Hao, Y., Pan, C., Sen, B., Novicki, M., Zhou, S., and Seftor, C. J.: Performance of the Ozone Mapping and Profiler Suite (OMPS) products, J. Geophys. Res.-Atmos., 119, 6181–6195,, 2014. a

Flynn, L. E., Zhang, Z., Mikles, V., Das, B., Niu, J., Beck, T. C., and Beach, E.: Algorithm Theoretical Basis Document for NOAA NDE OMPS Version 8 Total Column Ozone (V8TOz) Environmental Data Record (EDR) Version`1.0, ATBD, (last access: 1 October 2022), 2016. a

Fortuin, J. P. F. and Kelder, H.: An ozone climatology based on ozonesonde and satellite measurements, J. Geophys. Res.-Atmos., 103, 31709–31734,, 1998. a

Gelaro, R., McCarty, W., Suárez, M. J., Todling, R., Molod, A., Takacs, L., Randles, C. A., Darmenov, A., Bosilovich, M. G., Reichle, R., Wargan, K., Coy, L., Cullather, R., Draper, C., Akella, S., Buchard, V., Conaty, A., da Silva, A. M., Gu, W., Kim, G. K., Koster, R., Lucchesi, R., Merkova, D., Nielsen, J. E., Partyka, G., Pawson, S., Putman, W., Rienecker, M., Schubert, S. D., Sienkiewicz, M., and Zhao, B.: The modern-era retrospective analysis for research and applications, version 2 (MERRA-2), J. Climate, 30, 5419–5454,, 2017. a, b

Grainger, J. F. and Ring, J.: Anomalous Fraunhofer Line Profiles, Nature, 193, 762–762,, 1962. a

Herman, J., Huang, L., McPeters, R., Ziemke, J., Cede, A., and Blank, K.: Synoptic ozone, cloud reflectivity, and erythemal irradiance from sunrise to sunset for the whole earth as viewed by the DSCOVR spacecraft from the earth–sun Lagrange 1 orbit, Atmos. Meas. Tech., 11, 177–194,, 2018. a, b, c

Herman, J. R. and Celarier, E. A.: Earth surface reflectivity climatology at 340–380 nm from TOMS data, J. Geophys. Res.-Atmos., 102, 28003–28011,, 1997. a

Herman, J. R., Bhartia, P. K., Torres, O., Hsu, C., Seftor, C., and Celarier, E.: Global distribution of UV-absorbing aerosols from Nimbus 7/TOMS data, J. Geophys. Res.-Atmos., 102, 16911–16922,, 1997. a

Joiner, J. and Bhartia, P. K.: Accurate determination of total ozone using SBUV continuous spectral scan measurements, J. Geophys. Res.-Atmos., 102, 12957–12969,, 1997. a

Joiner, J. and Vasilkov, A. P.: First results from the OMI rotational Raman scattering cloud pressure algorithm, IEEE T. Geosci. Remote, 44, 1272–1282,, 2006. a

Joiner, J., Bhartia, P. K., Cebula, R. P., Hilsenrath, E., McPeters, R. D., and Park, H.: Rotational Raman scattering (Ring effect) in satellite backscatter ultraviolet measurements, Appl. Optics, 34, 4513,, 1995. a

King, M. D., Platnick, S., Menzel, W. P., Ackerman, S. A., and Hubanks, P. A.: Spatial and Temporal Distribution of Clouds Observed by MODIS Onboard the Terra and Aqua Satellites, IEEE T. Geosci. Remote, 51, 3826–3852,, 2013. a

Kleipool, Q. L., Dobber, M. R., de Haan, J. F., and Levelt, P. F.: Earth surface reflectance climatology from 3 years of OMI data, J. Geophys. Res.-Atmos., 113, D18308,, 2008. a, b

Koelemeijer, R. B. A.: A database of spectral surface reflectivity in the range 335–772 nm derived from 5.5 years of GOME observations, J. Geophys. Res.-Atmos., 108, 4070,, 2003. a

Koelemeijer, R. B. A. and Stammes, P.: Effects of clouds on ozone column retrieval from GOME UV measurements, J. Geophys. Res.-Atmos., 104, 8281–8294,, 1999. a

Labow, G. J., Ziemke, J. R., McPeters, R. D., Haffner, D. P., and Bhartia, P. K.: A total ozone-dependent ozone profile climatology based on ozonesondes and Aura MLS data, J. Geophys. Res.-Atmos., 120, 2537–2545,, 2015. a

Lamsal, L. N., Weber, M., Tellmann, S., and Burrows, J. P.: Ozone column classified climatology of ozone and temperature profiles based on ozonesonde and satellite data, J. Geophys. Res.-Atmos., 109, D20304,, 2004. a

Landgraf, J., Hasekamp, O., van Deelen, R., and Aben, I.: Rotational Raman scattering of polarized light in the Earth atmosphere: a vector radiative transfer model using the radiative transfer perturbation theory approach, J. Quant. Spectrosc. Ra., 87, 399–433,, 2004. a

Lerot, C., Van Roozendael, M., Lambert, J.-C., Granville, J., van Gent, J., Loyola, D., and Spurr, R. J. D.: The GODFIT algorithm: a direct fitting approach to improve the accuracy of total ozone measurements from GOME, Int. J. Remote Sens., 31, 543–550,, 2010. a

Lerot, C., Van Roozendael, M., Spurr, R. J. D., Loyola, D., Coldewey-Egbers, M. R., Kochenova, S., Van Gent, J., Koukouli, M., Balis, D., Lambert, J. C., Granville, J., and Zehner, C.: Homogenized total ozone data records from the European sensors GOME/ERS-2, SCIAMACHY/Envisat, and GOME-2/MetOp-A, J. Geophys. Res.-Atmos., 119, 1639–1662,, 2014. a, b

Levelt, P. F., Hilsenrath, E., Leppelmeier, G. W., van den Oord, G. H. J., Bhartia, P. K., Tamminen, J., de Haan, J. F., and Veefkind, J. P.: Science objectives of the ozone monitoring instrument, IEEE T. Geosci. Remote, 44, 1199–1208,, 2006. a

Liu, G., Tarasick, D. W., Fioletov, V. E., Sioris, C. E., and Rochon, Y. J.: Ozone correlation lengths and measurement uncertainties from analysis of historical ozonesonde data in North America and Europe, J. Geophys. Res.-Atmos., 114, D04112,, 2009. a

Loyola, D. G., Koukouli, M. E., Valks, P., Balis, D. S., Hao, N., Van Roozendael, M., Spurr, R. J. D., Zimmer, W., Kiemle, S., Lerot, C., and Lambert, J.-C.: The GOME-2 total column ozone product: Retrieval algorithm and ground-based validation, J. Geophys. Res.-Atmos., 116, D07302,, 2011. a

Lucht, W., Schaaf, C., and Strahler, A.: An algorithm for the retrieval of albedo from space using semiempirical BRDF models, IEEE T. Geosci. Remote, 38, 977–998,, 2000. a

Malicet, J., Daumont, D., Charbonnier, J., Parisse, C., Chakir, A., and Brion, J.: Ozone UV spectroscopy. II. Absorption cross-sections and temperature dependence, J. Atmos. Chem., 21, 263–273,, 1995. a

McPeters, R. D. and Labow, G. J.: Climatology 2011: An MLS and sonde derived ozone climatology for satellite retrieval algorithms, J. Geophys. Res.-Atmos., 117, D10303,, 2012. a

McPeters, R. D., Labow, G. J., and Logan, J. A.: Ozone climatological profiles for satellite retrieval algorithms, J. Geophys. Res.-Atmos., 112, D05308,, 2007. a, b, c

Munro, R., Lang, R., Klaes, D., Poli, G., Retscher, C., Lindstrot, R., Huckle, R., Lacan, A., Grzegorski, M., Holdak, A., Kokhanovsky, A., Livschitz, J., and Eisinger, M.: The GOME-2 instrument on the Metop series of satellites: instrument design, calibration, and level 1 data processing – an overview, Atmos. Meas. Tech., 9, 1279–1301,, 2016. a

NASA/LARC/SD/ASDC: DSCOVR EPIC Level 2 O3SO2AI, NASA Langley Atmospheric Science Data Center DAAC [data set],, 2018. a

Platt, U.: Air Monitoring by Differential Optical Absorption Spectroscopy, in: Encyclopedia of Analytical Chemistry, iii, John Wiley & Sons, Ltd, Chichester, UK, 1–28,, 2017. a

Rodgers, C. D.: Inverse Methods for Atmospheric Sounding, Series on Atmospheric, Oceanic and Planetary Physics, Vol. 2, edited by: Taylor, F. W., World Scientific,, 2000. a

Schaaf, C. B., Liu, J., Gao, F., and Strahler, A. H.: Aqua and Terra MODIS Albedo and Reflectance Anisotropy Products, in: Land Remote Sensing and Global Environmental Change, edited by: Ramachandran, B., Justice, C., and Abrams, M., Remote Sensing and Digital Image Processing, Vol. 11, Springer, New York, NY,, 2010. a

Sofieva, V. F., Tamminen, J., Kyrölä, E., Mielonen, T., Veefkind, P., Hassler, B., and Bodeker, G. E.: A novel tropopause-related climatology of ozone profiles, Atmos. Chem. Phys., 14, 283–299,, 2014. a

Spurr, R. J.: VLIDORT: A linearized pseudo-spherical vector discrete ordinate radiative transfer code for forward model and retrieval studies in multilayer multiple scattering media, J. Quant. Spectrosc. Ra., 102, 316–342,, 2006. a, b, c

Spurr, R. J. D., de Haan, J., van Oss, R., and Vasilkov, A. P.: Discrete-ordinate radiative transfer in a stratified medium with first-order rotational Raman scattering, J. Quant. Spectrosc. Ra., 109, 404–425,, 2008. a

Stajner, I., Wargan, K., Pawson, S., Hayashi, H., Chang, L.-P., Hudman, R. C., Froidevaux, L., Livesey, N., Levelt, P. F., Thompson, A. M., Tarasick, D. W., Stübi, R., Andersen, S. B., Yela, M., König-Langlo, G., Schmidlin, F. J., and Witte, J. C.: Assimilated ozone from EOS-Aura: Evaluation of the tropopause region and tropospheric columns, J. Geophys. Res.-Atmos., 113, D16S32,, 2008. a

Tilstra, L. G., Tuinder, O. N. E., Wang, P., and Stammes, P.: Surface reflectivity climatologies from UV to NIR determined from Earth observations by GOME-2 and SCIAMACHY, J. Geophys. Res.-Atmos., 122, 4084–4111,, 2017. a, b, c

Torres, O., Bhartia, P. K., Herman, J. R., Ahmad, Z., and Gleason, J.: Derivation of aerosol properties from satellite measurements of backscattered ultraviolet radiation: Theoretical basis, J. Geophys. Res.-Atmos., 103, 17099–17110,, 1998. a

Torres, O., Tanskanen, A., Veihelmann, B., Ahn, C., Braak, R., Bhartia, P. K., Veefkind, P., and Levelt, P.: Aerosols and surface UV products from Ozone Monitoring Instrument observations: An overview, J. Geophys. Res.-Atmos., 112, D24S47,, 2007. a, b

Van Roozendael, M., Loyola, D., Spurr, R., Balis, D., Lambert, J.-C., Livschitz, Y., Valks, P., Ruppert, T., Kenter, P., Fayt, C., and Zehner, C.: Ten years of GOME/ERS-2 total ozone data – The new GOME data processor (GDP) version 4: 1. Algorithm description, J. Geophys. Res.-Atmos., 111, D14311,, 2006. a

Van Roozendael, M., Spurr, R. J. D., Loyola, D., Lerot, C., Balis, D., Lambert, J. C., Zimmer, W., Van Gent, J., Van Geffen, J., Koukouli, M., Granville, J., Doicu, A., Fayt, C., and Zehner, C.: Sixteen years of GOME/ERS-2 total ozone data: The new direct-fitting GOME Data Processor (GDP) version 5 Algorithm description, J. Geophys. Res.-Atmos., 117, D03305,, 2012. a

Vasilkov, A. P., Joiner, J., Spurr, R. J. D., Bhartia, P. K., Levelt, P., and Stephens, G.: Evaluation of the OMI cloud pressures derived from rotational Raman scattering by comparisons with other satellite data and radiative transfer simulations, J. Geophys. Res.-Atmos., 113, D15,, 2008. a

Veefkind, J., de Haan, J., Brinksma, E., Kroon, M., and Levelt, P.: Total ozone from the ozone monitoring instrument (OMI) using the DOAS technique, IEEE T. Geosci. Remote, 44, 1239–1244,, 2006. a

Veefkind, J. P., Aben, I., McMullan, K., Förster, H., de Vries, J., Otter, G., Claas, J., Eskes, H. J., de Haan, J. F., Kleipool, Q., van Weele, M., Hasekamp, O., Hoogeveen, R., Landgraf, J., Snel, R., Tol, P., Ingmann, P., Voors, R., Kruizinga, B., Vink, R., Visser, H., and Levelt, P. F.: TROPOMI on the ESA Sentinel-5 Precursor: A GMES mission for global observations of the atmospheric composition for climate, air quality and ozone layer applications, Remote Sens. Environ., 120, 70–83,, 2012. a

Vountas, M., Rozanov, V., and Burrows, J.: Ring Effect: Impact Of Rotational Raman Scattering On Radiative Transfer In Earth's Atmosphere, J. Quant. Spectrosc. Ra., 60, 943–961,, 1998. a

Wagner, T., Beirle, S., Deutschmann, T., and Penning de Vries, M.: A sensitivity analysis of Ring effect to aerosol properties and comparison to satellite observations, Atmos. Meas. Tech., 3, 1723–1751,, 2010. a

Wargan, K., Pawson, S., Olsen, M. A., Witte, J. C., Douglass, A. R., Ziemke, J. R., Strahan, S. E., and Nielsen, J. E.: The global structure of upper troposphere-lower stratosphere ozone in GEOS-5: A multiyear assimilation of EOS Aura data, J. Geophys. Res.-Atmos., 120, 2013–2036,, 2015. a

Wargan, K., Labow, G., Frith, S., Pawson, S., Livesey, N., and Partyka, G.: Evaluation of the Ozone Fields in NASA's MERRA-2 Reanalysis, J. Climate, 30, 2961–2988,, 2017. a

Wassmann, A., Borsdorff, T., aan de Brugh, J. M. J., Hasekamp, O. P., Aben, I., and Landgraf, J.: The direct fitting approach for total ozone column retrievals: a sensitivity study on GOME-2/MetOp-A measurements, Atmos. Meas. Tech., 8, 4429–4451,, 2015. a

Wei, J. C., Pan, L. L., Maddy, E., Pittman, J. V., Divarkarla, M., Xiong, X., and Barnet, C.: Ozone profile retrieval from an advanced infrared sounder: experiments with tropopause-based climatology and optimal estimation approach, J. Atmos. Ocean. Tech., 27, 1123–1139,, 2010. a

Wellemeyer, C. G., Taylor, S. L., Seftor, C. J., Mcpeters, R. D., and Bhartia, P. K.: A correction for total ozone mapping spectrometer profile shape errors at high latitude, J. Geophys. Res.-Atmos., 102, 9029,, 1997. a, b

Yang, K. and Liu, X.: Ozone profile climatology for remote sensing retrieval algorithms, Atmos. Meas. Tech., 12, 4745–4778,, 2019. a, b, c, d, e, f, g

Yang, K., Fleig, A., Wolfe, R., and Nishihama, M.: MODIS band-to-band registration, in: IGARSS 2000. IEEE 2000 International Geoscience and Remote Sensing Symposium. Taking the Pulse of the Planet: The Role of Remote Sensing in Managing the Environment. Proceedings (Cat. No.00CH37120), IEEE 2, Honolulu, HI, USA, 24–28 July 2000, 887–889,, 2000. a

Yang, K., Bhartia, P. K., Wellemeyer, C. G., Qin, W., Spurr, R. J. D., Veefkind, J. P., and de Haan, J. F.: Application of spectral fitting method to GOME and comparison with OMI‐DOAS and TOMS‐V8 total ozone, in: Proceedings of the XX Quadrennial Ozone Symposium, edited by Zerefos, C. S., International Ozone Commission, 1–8 June 2004, Kos, Greece, 510–511, (last access: 1 October 2022), 2004. a

Yang, K., Krotkov, N. A., Krueger, A. J., Carn, S. A., Bhartia, P. K., and Levelt, P. F.: Retrieval of large volcanic SO2 columns from the Aura Ozone Monitoring Instrument: Comparison and limitations, J. Geophys. Res.-Atmos., 112, D24S43,, 2007. a

Yang, K., Krotkov, N. A., Krueger, A. J., Carn, S. A., Bhartia, P. K., and Levelt, P. F.: Improving retrieval of volcanic sulfur dioxide from backscattered UV satellite observations, Geophys. Res. Lett., 36, L03102–L03102,, 2009a. a

Yang, K., Liu, X., Krotkov, N. a., Krueger, A. J., and Carn, S. a.: Estimating the altitude of volcanic sulfur dioxide plumes from space borne hyper-spectral UV measurements, Geophys. Res. Lett., 36, L10803–L10803,, 2009b. a, b

Yang, K., Liu, X., Bhartia, P. K., Krotkov, N. A., Carn, S. A., Hughes, E. J., Krueger, A. J., Spurr, R. J. D., and Trahan, S. G.: Direct retrieval of sulfur dioxide amount and altitude from spaceborne hyperspectral UV measurements: Theory and application, J. Geophys. Res.-Atmos., 115, D00L09,, 2010. a, b, c, d

Yang, K., Dickerson, R. R., Carn, S. A., Ge, C., and Wang, J.: First observations of SO2 from the satellite Suomi NPP OMPS: Widespread air pollution events over China, Geophys. Res. Lett., 40, 4957–4962,, 2013. a

Yang, K., Carn, S. a., Ge, C., Wang, J., and Dickerson, R. R.: Advancing Measurements of Tropospheric NO2 from Space: New Algorithm and First Global Results from OMPS, Geophys. Res. Lett., 41, 4777–4786,, 2014.  a

Yang, Y., Meyer, K., Wind, G., Zhou, Y., Marshak, A., Platnick, S., Min, Q., Davis, A. B., Joiner, J., Vasilkov, A., Duda, D., and Su, W.: Cloud products from the Earth Polychromatic Imaging Camera (EPIC): algorithms and initial evaluation, Atmos. Meas. Tech., 12, 2019–2031,, 2019. a

Short summary
This paper describes the algorithm for O3 and SO2 retrievals from DSCOVR EPIC. Algorithm advances, including the improved O3 profile representation and the regulated direct fitting inversion technique, improve the accuracy of O3 and SO2 from the multi-channel measurements of DSCOVR EPIC. A thorough error analysis is provided to quantify O3 and SO2 retrieval uncertainties due to various error sources and simplified algorithm physics treatments.