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

. On board the Deep Space Climate Observatory (DSCOVR), the ﬁrst Earth-observing satellite at the L1 point (the ﬁrst 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 reﬂected radiances in 10 discrete spectral channels, four of which are in the ultraviolet (UV) range. These UV bands are selected primarily for total ozone (O 3 ) and aerosol retrievals based on heritage algorithms developed for the series of Total Ozone Mapping Spectrome-ters (TOMS). These UV measurements also provide sensitive detection of sulfur dioxide (SO 2 ) 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 ﬁtting (DVCF) algorithm used for retrieving total vertical columns of O 3 and SO 2 from DSCOVR EPIC. This paper describes algorithm advances, including an improved O 3 proﬁle representation that enables proﬁle ad-justments from multiple spectral measurements and the spatial optimal estimation (SOE) scheme that reduces O 3 artifacts resulting from EPIC’s band-to-band misregistrations. Furthermore, this paper discusses detailed error analyses and presents intercomparisons with correlative data to validate O 3 and SO 2 retrievals from EPIC.


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×10 6 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) chargedcoupled 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.
Published by Copernicus Publications on behalf of the European Geosciences Union. 5878 X. Huang and K. Yang: O 3 and SO 2 ATBD 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 km 2 at the image center of the sunlit disk. The effective footprint size is about 10 × 10 km 2 , 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 km 2 , 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 km 2 ) on a series of satellites, the Scanning Imaging Absorption Spectrometer for Atmospheric Cartography (SCIAMACHY; 60 × 30 km 2 ; Bovensmann et al., 1999) on ESA's ENVIronmental SATellite (ENVISAT), the Ozone Mapping and Profiler Suite Nadir Mapper (OMPS-NM; 50 × 50 km 2 ; 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 km 2 ), Metop-B (80 × 40 km 2 ), and Metop-C (80 × 40 km 2 ). Though it is slightly larger than the nadir footprint of the Ozone Monitoring Instrument (OMI; 13 × 24 km 2 ; Levelt et al., 2006) on Aura and the OMPS-NM (17 × 13 km 2 ; Flynn et al., 2016) on NOAA-20, as well as much bigger than that of the TROPOspheric Monitoring Instrument (TROPOMI; 5.5×3.5 km 2 ; Veefkind et al., 2012) on the ESA Sentinel-5 Precursor (S5P), EPIC's spatial resolution is sufficiently high to map small-scale O 3 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).
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 (O 3 ) retrievals. These UV bands also provide sensitive detection of sulfur dioxide (SO 2 ) 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 O 3 and SO 2 retrievals achieved by applying the DVCF algorithm to spectral UV radiance measurements of DSCOVR EPIC. Lastly, we validate the DSCOVR EPIC O 3 and SO 2 through intercomparisons with correlative data.

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 L m (in units of W · sr −1 · m −2 · nm −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 where F (λ) (in units of W · m −2 · nm −1 ) is the monochromatic spectral solar irradiance, and I TOA (λ) is the Sunnormalized monochromatic top-of-the-atmosphere (TOA) radiance (in units of sr −1 ) for a wavelength λ (in units of nm). The Sun-normalized measured radiance I M for a spectral band is defined as I M = L M /F M , where F M = 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 I M , 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, I TOA (λ) and other spectral-dependent quantities are hereafter used to denote flux-weighted bandpass averages, with λ representing the characterized wavelength of the spectral band.
To reach a sensor at TOA, sunlight photons are either backscattered 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., O 3 and SO 2 ) 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 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. Figure 2. Filter 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). 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 altituderesolved air mass factor (AMF, m z ). 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, I TOA , 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 λ, I TOA can be expressed as the sum of two contributions, where I a consists of solar photons scattered once or more by molecules and particles in the atmosphere without interacting with the underlying surface, and I s represents solar photons reflected at least once or multiple times by the underlying surface.

Path radiance
I a 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) (Dave, 1964) in the direction specified by the view zenith angle (θ v ), attenuated (e −t/µ , where µ = cos θ v ) 5880 X. Huang and K. Yang: O 3 and SO 2 ATBD 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: 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 I a 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 I a 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, Dave, 1964, or VLIDORT, Spurr, 2006 as a function of wavelength for a molecular atmosphere with the O 3 profile X 1 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 O 3 ). O 3 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 σ (O 3 ) 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 O 3 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 O 3 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 (ψ = /I a ) of EPIC bands 1 and 2 for two different observation geometries and a midlatitude O 3 profile labeled as X 1 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.
The measurement sensitivity to a thin molecular absorber layer is equal to the product of the absorption cross-sections (σ ) and the mean path length (m a ) of photons passing through the layer, where m a = −∂ ln I a /∂τ 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, m G = 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 m a decreases rapidly as the layer descends nearing the surface due to fewer photons reaching the lower atmosphere, while m a approaches m G 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), m a of the low zenith geometry usually exceeds m G due to a significant fraction of photons undergoing multiple scattering below and within UTLS, while m a 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 m a is shorter for a wavelength with stronger O 3 absorption, which reduces the number of photons reaching the lower atmosphere. The variation of m a with a changing altitude signifies the path radiance dependence on the absorber profile. The path radiance fractional change due to profile change ( X = X 2 − X 1 ) can be expressed as where T z is the atmospheric temperature, and X 1 (z) and X 2 (z) are absorber concentration at altitude z. Figure 3b illustrates the change in path radiance caused by an O 3 profile change while keeping its total vertical column the same: lowering the O 3 profile (e.g., X 1 to X 2 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 O 3 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 Figure 3. Sample results from RT simulations for a molecular atmosphere with O 3 profiles X 1 and X 2 in panel (a). Both X 1 and X 2 are midlatitude zone (30 • ≤ latitude ≤ 60 • ) climatological O 3 profiles with the same total vertical column of 275 Dobson units, where 1 DU = 2.69 × 10 16 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 I a (X 1 ) for the low and high zenith geometries and their fractional changes ( I a /I a ) when the O 3 profile is changed to X 2 . (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 (m a ) 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, m G . 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 O 3 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 O 3 measurement.

Surface reflection
The path radiance I a includes backscattered photons that are independent of the underlying surface, while the surface contribution to TOA radiance, I s (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 I s can be accurately predicted with RT modeling. For a Lambertian surface, which reflects radiation isotropically independent of the incident direction, the surface radiance I s can be expressed as (Dave, 1964) where r s 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 S b 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, I s , may be described as the once-reflected radiance (T ↓ r s T ↑ ), 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 − r s S b ).
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 Bandeen, 1970) 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 r s that reproduces the radiance I s 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.
Reflection of UV sunlight from natural surfaces has long been measured by instruments on board satellites in Sunsynchronous 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 r s from a measured radiance I M , the atmospheric path radiance I a , transmissions T ↓ and T ↑ , and reflectance S b for a spectral band are calculated for a molec- where I s = I M − I a . 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 , GOME-1 at 335-772 nm from 1995-2000 (Koelemeijer, 2003), 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 cloudand 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 Reynolds, 1971;Green, 1980, 1981;Feister and Grewe, 1995), 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 O 3 absorption) wavelengths to estimate the GLERs at short (strong O 3 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 I s 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, m s = −∂ ln I s /∂τ 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 O 3 (∼ 90 %) is located in the stratosphere, the Lambertian treatment does not introduce a significant AMF error in total O 3 absorption.
As described above, UV reflectivities for most natural surfaces are quite low (GLER < 0.1); therefore, the surface contributions I s are typically much smaller than (< 10 % at 317.5 nm) the path radiance I a (see Fig. 6). In modeling a measured radiance I M , an error in surface radiance I s is compensated for with the path radiance I a . The uncertainty of extrapolated GLER is usually less than 1 %, corresponding to a less than 1 % error in I s and hence less than 0.1 % error in the path radiance I a . Furthermore, the AMF error due to the Lambertian treatment of an anisotropic surface is insignificant, since the combined mean photon path lengths, m z = −∂ ln I TOA /∂τ z = (I a m a + I s m s )/I TOA , contain minor contributions from surface radiance I s . 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 I s (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 O 3 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 O 3 absorption is suppressed for the low SZA and VZA observations.

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 sur-face 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 the weighted sum of two independent contributions I g and I c . Here I g is the radiance from the cloud-free portion of the IFOV containing a Lambertian surface of reflectivity R g at pressure p g . Similarly, I c is from the cloudy portion, and f c is the cloud fraction and R c the reflectivity of the Lambertian surface at pressure p c . The MLER model can reproduce measured radiances I m through the determination of cloud fraction f c . First, the scene LER r s at surface pressure p g is estimated using Eq. (6). If r s is less than or equal to the climatological LER value R g (e.g., Kleipool et al., 2008), this IFOV is treated as a particle-free scene (f c = 0). If r s is greater than or equal to the LER value for cloud R c = 0.8 (Koelemeijer and Stammes, 1999;Ahmad et al., 2004), this IFOV is treated as fully cloud-covered (f c = 1). When r s is between R g and R c , the cloud fraction is inverted from Eq. (8), which yields In the case of f c = 0 or 1, surface LER r g or cloud LER r c is determined using Eq. (6) to ensure that modeled radiance I TOA is equal to the measurement I M . Figure 8 shows cloud fractions (f c ) as a function of wavelength for several examples of particle-laden atmospheres. The radiance intensity scattered from atmospheric particles varies smoothly with wavelength without highfrequency spectral structures. For instance, the contributions to TOA radiances (I TOA ) from backscattering by meteorological clouds change smoothly and slowly with wavelength (see Fig. 8, the CLD curve). The selection of R c = 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 f c has a small spectral variation (i.e., f c nearly the same for different wavelengths) for most cloudy observations. The small and smooth change of f c with wavelengths allows its extrapolation to provide a reliable estimate of f c 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 I TOA due to the addition of aerosols and hence the cloud fraction (f c or the surface LER, r g ) are smooth in wavelength (see, e.g., Fig. 8, with smooth curves for weakly absorbing sulfatebased aerosols -SLFs, carbonaceous aerosols from biomass Figure 8. Four examples of cloud fractions (f c ) 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; Deirmendjian, 1969) 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, R g , p g , R c , and p c , as well as the angles (θ s , θ v , and φ) that specify the viewing illumination geometry.
burning -BIO, and mineral dust -DST; Torres et al., 2007). Therefore, f c (when f c > 0, from Eq. 9) or r g (when f c = 0, from Eq. 6) determined at longer wavelengths where atmospheric absorption is weak may be linearly extrapolated to O 3 -sensitive wavelengths for estimation of contributions to TOA radiance from surface reflection and particle backscattering (referred to as the r g f c extrapolation method hereafter).
The UV aerosol index (AI; 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 c l used in the r g f c extrapolation scheme. Algebraically, AI is calculated as the N-value (defined as −100log 10 I ) difference between the modeled (I TOA ) and measured (I M ) radiances at a wavelength λ.
= 100 c l λ ∂log 10 I TOA (λ, R) ∂R R=R e Here, the modeled radiance I TOA (λ, R e ) is calculated for a molecular atmosphere with an estimated reflectivity parameter R e , which may be the LER value r s or the MLER cloud fraction f c 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 I M (λ) = I TOA (λ, R m ) = I TOA (λ, R e + R), since the reflectivity parameter R m is derived from I M (λ) and R = R m − R e = c l λ, we arrive at Eq. (11) from the definition of AI in Eq. (10). In short, the spectral slope c l 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 I TOA is modeled for a Rayleighscattering-only atmosphere over a Lambertian surface. To capture the spectral slope of the r s f c extrapolation scheme, we switch the LER treatment with the MLER modeling of I TOA for AI calculation. The resulting MLER AI is usually higher than the corresponding LER AI when f c > 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 (r g or f c ) follows the linear relationship among different wavelengths. In reality, the spectral dependence of natural surface reflection (r g ) or particle scattering and absorption (f c ) is nonlinear, though moderately as exemplified in Figs. 4b and 8. Therefore, r g f c extrapolation yields small errors in r g or f c at the extrapolated wavelengths. The radiance uncertainties associated with the r g f c 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 O 3 absorption, which occurs mostly in the stratosphere. This is similar to how the Lambertian treatment of surface reflection works for the estimation of total O 3 absorption (see Sect. 2.2).
The MLER treatment relies on a few adjustable parameters, including the cloud fraction f c and cloud pressure p c , to model a vast range of conditions encountered in remote sensing of Earth's atmosphere. The cloud fraction f c , obtained directly from radiance measurements using Eq. (9), provides an estimate of the cloud amount in an IFOV. The pressure p c 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 p c located too high or too low from the optical centroid pressure (OCP; Joiner and Vasilkov, 2006;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 O 2 A band  are usually located within the particle vertical distributions and therefore used to set the cloud pressures p c for processing EPIC observations.
The use of OCP for p c 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 p c 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 O 3 , such as NO 2 and SO 2 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.

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 Spurr, 1997;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 Ring, 1962).
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 = (I RRS − I ELA )/I ELA , where I ELA is the TOA radiance calculated assuming all molecular scattering is elastic, while I RRS 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 . Since RRS is weakly dependent on polarization, a scalar radiative transfer model, from which both I ELA and I RRS are calculated without including radiation polarization, can accu- × 100 for EPIC UV bands as a function zenith angle at two surface pressures p g = 0.7 ATM and p g = 1.0 ATM, with a surface albedo of r g = 0.1. rately provide filling-in factors (Landgraf et al., 2004;Wagner et al., 2010).
The filling-in factors provide estimates of the modeling errors in I TOA 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 %) O 3 retrieval accuracy. The filling-in factors ( ), modeled using a scaler code (like LIDORT-RRS), may be used to correct the results (I TOA ) from vector radiative transfer codes (e.g., Dave, 1964;Spurr, 2006) that perform elastic modeling only, i.e., the RRS-corrected TOA radiance = I TOA ( + 1).

O and temperature vertical profiles
As shown in Sect. 2.1, the O 3 vertical distribution or profile directly affects the magnitude of a measured radiance in the spectral region with significant O 3 absorption. Hence, the interpretation of radiance change due to O 3 absorption requires some knowledge of its profile. In general, the retrieval of quantitative information about a gaseous absorber (such as O 3 and SO 2 ) 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 O 3 profile model for remote sensing retrieval algorithms and its improvements over the model commonly used by other total O 3 algorithms.
Ozone is naturally present throughout the atmosphere, and its spatial and temporal distribution is controlled by atmospheric processes of O 3 production, destruction, and transport. The O 3 distribution exhibits a high abundance of O 3 in the stratosphere and a minor portion (∼ 10 %) in the tropo-sphere, with the peak O 3 concentration occurring at a lower altitude as the latitude increases towards the poles. These characteristics are well-captured by O 3 profile climatologies (e.g., Fortuin and Kelder, 1998;McPeters et al., 2007;McPeters and Labow, 2012), which provide the mean and variance of O 3 vertical distribution as a function of latitude and calendar month. These climatologies also reveal that the O 3 profile has the highest variability in the upper troposphere and lower stratosphere (UTLS), contributing the most to the natural variations in total O 3 . This high O 3 variability is the consequence of atmospheric movements that blend air masses with different O 3 concentrations, such as uplifting of O 3 -poor air in the troposphere or lowering of O 3 -rich air in the stratosphere resulting from the rise and fall of the tropopause. Predictors of O 3 profile shape, including tropopause pressure and total O 3 columns, are developed to capture the dynamical influences on O 3 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 Wellemeyer, 2002;Lamsal et al., 2004;Labow et al., 2015) O 3 profile climatologies.
To improve the representation of the O 3 profile, we construct both tropopause-dependent and total-columndependent climatologies using the Modern-Era Retrospective Analysis for Research and Applications Version 2 (MERRA-2; Bosilovich et al., 2015;Gelaro et al., 2017) O 3 record between 2005 and 2016. The total-column-dependent climatology, named M2TCO3, is more appropriate for use as the O 3 profile model needed by a total O 3 algorithm, as it is generally more reliable than the tropopause-dependent version in prescribing realistic O 3 profiles (Yang and Liu, 2019). Figure 11 compares daytime M2TCO3 (Yang and Liu, 2019; referred to as M2TCO3 hereafter) and TOMS-V8 profiles for 2 months and four latitude zones, illustrating the similarities and differences between the two O 3 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. Fig. 11), a strong dependence on latitude exhibiting lower altitudes of O 3 concentration peaks at higher latitudes for similar total columns, and characteristic dependence on the total column, which gets smaller with a higher O 3 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 Wellemeyer, 2002), thus ignoring the significant seasonal dependence. An additional deficiency of TOMS-V8 contributing to the differences is its inadequate representation of latitude-dependent O 3 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 O 3 profile representation.

September in
In short, M2TCO3 better captures the dynamical changes and spatiotemporal variations in O 3 profiles with higher resolutions in total O 3 column (25 DU), latitude (10 • ), and time (monthly). Taking into account the substantial change in atmospheric O 3 over a long time, M2TCO3 is more accurate in representing atmospheric O 3 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 O 3 profile model for total O 3 retrieval from EPIC.
The M2TCO3 climatology contains not only mean profiles that represent the likely O 3 vertical distributions, but also the modal O 3 adjustment profiles that specify the probable deviations from the means. These modal profiles are determined from the O 3 profile covariance statistics, as illustrated in Fig. 12, showing examples of M2TCO3 climatological O 3 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 O 3 profile X is expressed as where X m (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. e k (v) is the kth modal profile, γ k the kth coefficient, and p the number of e k (v), with a maximum equal to the number of levels used to represent an O 3 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 O 3 profile X, which deviates invariably from the mean X m , can be accurately represented using Eq. (12) with a small number of expansion coefficients γ k . Much like the mean the profile X m represents the most probable vertical distribution of O 3 : the modal profiles, {e k , 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 O 3 profile information contained in multi-spectral measurements to improve the O 3 profile representation by determining one or more linear expansion coefficients {γ k , k = 1. . .}. Note that for most total ozone algorithms, the O 3 profile representation is limited to the climatological mean only, equivalent to restricting γ k = 0 for all k in Eq. (12). The total column is a good predictor of an O 3 profile that is especially accurate for the shape in the stratosphere but less so in the troposphere. Tropospheric O 3 exhibits characteristic spatiotemporal distribution, which is captured in the MERRA-2 tropospheric O 3 climatology (Yang and Liu, 2019). To better represent the O 3 profile, the tropospheric part of a column-dependent M2TCO3 profile, X m , 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 X m in Eq. (12) has its tropospheric part tied to the spatiotemporally varying climatological tropospheric column, to which the tropospheric column of the mean X m profile (obtained by averaging over the different column amounts) is matched.
In addition to knowledge of profiles of light-absorbing trace gases, such as O 3 and SO 2 , radiative transfer modeling of measured radiance requires knowledge of the atmospheric temperature profile because the absorption crosssections of these trace gases significantly depend on temperature. For total O 3 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 Liu, 2019). This temperature climatology provides mean temperature profiles corresponding to the climatological O 3 profiles, capturing the dependence of the temperature profile on season and location, as well as the variation of temperature with the O 3 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 O 3 profiles.

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 instrumen-tal 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 O 3 and SO 2 , 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 O 3 retrieval, Eq. (12) embodies the profile constraint, while the MLER model with r s f c 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 O 3 retrieval, consist of total O 3 column , a number (p) of modal expansion coefficients {γ k , k = 1. . .p}, surface LER (r s ) or cloud fraction (f c ), and a number (q) of polynomial coefficients {c l , l = 1. . .q} of the r s f c 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 O 3 profile, and q = 0 for the spectral invariant reflectivity parameter.

Exact solution
Conceptually, the inversion is to find the state vector (x) that satisfies a set of m simultaneous equations, { y i = 0, i = 1. . .m}, one for each spectral band difference, y i = ln I M (λ i ) − ln I TOA (x, λ i ), between the radiance measurement I M and the forward modeling I TOA : here λ i the wavelength that characterizes the ith (1 ≤ i ≤ m) spectral band and y i the residual of this band. In matrix form, the m simultaneous equations can be expressed as where y is residual column vector { y i , i = 1. . .m}. Since the forward mapping I TOA (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 x L to linearize the equation between residuals and the state vector.
where x j and x Lj are the j th components of x and x L , respectively, x j = x j − x Lj the j th components of state adjustment vector, and K ij = ∂ ln I TOA (x,λ i ) ∂x j | x=x L the Jacobian, also known as the weighting function for the retrieval parameter x j at spectral band λ i . The m residual elements, each written in Eq. (14), can be expressed in matrix form as where y L is the column vector {ln I M (λ i ) − ln I TOA (x L , λ i ), i = 1. . .m}, x = x − x L the state adjustment vector, and K the m × n Jacobian matrix with the {K ij , i = 1. . .m, j = 1. . .n} as its elements. Putting Eq. (15) into Eq. (13) yields 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 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 O 3 algorithm (Bhartia and Wellemeyer, 2002). The TOMS-V8 algorithm determines the two-component state vector, x = { , r s or f c }, from radiance measurements of two spectral bands: one with low O 3 sensitivity to estimate the MLER parameter (r s or f c ) and the other with high O 3 sensitivity to derive total O 3 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 O 3 and SO 2 columns from EPIC UV measurements.

Direct fitting
Since spectral measurements have errors and m = n 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 µ 2 i . 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 x L . Substituting y (Eq. 15) into Eq. (18), we minimize this cost function to obtain the state adjustment vector which is the solution of linear weighted least-square regression. Here, G DF = (K T S −1 K) −1 K T S −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 O 3 vertical column (Joiner and Bhartia, 1997;Yang et al., 2004;Lerot et al., 2014), the combination of total O 3 and SO 2 vertical columns (Yang et al., 2007(Yang et al., , 2009a(Yang et al., , 2013, the combination of O 3 and altitude-resolved SO 2 vertical columns (Yang et al., 2009b(Yang et al., , 2010, and stratospheric and tropospheric NO 2 vertical columns (Yang et al., 2014). This algorithm is named DVCF to contrast with the DOAS (differential optical absorption spectroscopy) method (Platt, 2017), 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 O 3 column ( ) and an expansion coefficient (γ k ) of differential profile e k (see Eq. 12), or the SO 2 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., 2009bYang et al., , 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.

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 where x a is the a priori state vector and S a the a priori state vector covariance matrix. The first term on the righthand 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 x a . 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 (Rodgers, 2000) to minimize the cost function, Eq. (21) yields the a posterior state adjustment vector at the Lth iteration: where x aL = x a − x L . Inserting I n = (K T S −1 K)(K T S −1 K) −1 , an n × n identity matrix, in front of the term K T S −1 y L 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 x 0 , and a frequent selection is the a priori state vector: x 0 = x a . 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 G DF y L (see Eq. 20, which is derived without any a priori constraint) and the difference between the a priori and the previous state vectors x aL weighted by matrices K T S −1 K and S −1 a , respectively. For the state vector component with a strong a priori constraint, i.e., a small variance in S a , 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 S a , 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 (S −1 a + K T S −1 K) −1 (see Eq. 23) and is thus less than or equal to the corresponding a priori variance in the a priori S a 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 S a . To further reduce the dependence on the a priori state vector, it is updated at each iteration with the linearization point, setting x a = x L , and hence Eq. (22) becomes where G = S −1 a + K T S −1 K −1 K T S −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 S a such that S −1 a → 0, G becomes G DF 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 (S a ) 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 O 3 and SO 2 retrievals.

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 O 3 and SO 2 absorptions (see Fig. 13), containing information that allows the retrieval of total O 3 and SO 2 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 O 3 and SO 2 absorption such that changes in O 3 and SO 2 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 (O 3 sensitive) channels, accomplished with the r s f c extrapolation scheme (see Sect. 2.3). The reflectivity spectral slope c l of this extrapolation is proportional to the AI (see Eq. 11). The reflectivity parameter (R) is either the LER value r s estimated from Eq. (6) or the MLER cloud fraction f c from Eq. (9) depending on the value of f c (R = r s when f c = 0, R = f c when f c > 0), and its spectral slope is calculated as c l = (R 4 − R 3 )/(λ 4 − λ 3 ).
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.

Reflectivity correction by spatial optimal estimation (SOE)
The estimation of the O 3 column from EPIC radiance measurements requires accurate reflectivity information on the underlying surface, which is extrapolated from the reflec- tivity 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 selfrotation and spacecraft jitter, different spectral images record slightly different (i.e., rotated) sunlit hemispheres. The geolocation procedure of EPIC (Blank, 2019) 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, r s f c 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 O 3 -sensitive band is to derive it from the radiance measurement of this band with an optimally estimated total O 3 column from the nearby O 3 distribution. This O 3 estimation is attainable because an actual spatial distribution of the total O 3 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 (R B ) at EPIC band B by minimizing the cost function that embodies the a priori knowledge of R B and O 3 spatial distribution. The first part of cost function supports a smooth (i.e., homogeneous) O 3 distribution, while the sec-ond part penalizes the difference between R B and its a priori value, which is the extrapolated reflectivity (R E ) from the longer wavelength EPIC bands. Hence the cost function is written as subject to the measurement constraints {I M (λ B , i) = I TOA ( (i), R B (i), λ B ), and IOFV index i = 1. . .n, the size of the IFOV group}, which is linearized to become 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 i-j pair. is the average O 3 column for the group. Given R E , which is the band-B reflectivity extrapolated from the longer wavelength bands, the total O 3 column (R E ) is retrieved from band-B radiance measurement using the exact solution method (see Sect. 4.1), and the associated O 3 profile is the column-dependent M2TCO3 climatological profile X m ( ). The equation of measurement constraint (Eq. 28) describes a positive (since S > 0 usually) linear relationship between total O 3 column and the surface reflectivity (in the neighborhood of R E ), increasing R B requires more O 3 absorption to maintain I M = I TOA .
Minimizing only the first rhs term of Eq. (26) leads to the same O 3 column for all the IFOVs (i.e., { (i) = , i = 1. . .n}), while minimizing only the second term makes R B = R E for each IFOV. The constants α and β are weights to respectively emphasize the smoothness of O 3 spatial distribution and the closeness of reflectivity between extrapolation and estimation. In the SOE scheme, weights are α = 0 and β = 1 for the traditional O 3 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 R B (i) to minimize each component ϒ i . The solution R B (i) that minimizes ϒ i is found by solving this equation: which yields From Eqs. (29) to (30), only the n nearby IFOVs are included, i.e., wt(i, j ) = 1 for i-j separation within a few (< 4) adjacent IFOVs; otherwise, wt(i, j ) = 0. At the start of iteration, { (j ) = (j, R E ), 1. . .n}, and they are then updated using Eq. (28) with R B (i) from Eq. (30) for the next iteration, which stops until changes in {R B (i), i = 1. . .n} 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 O 3 map (Fig. 14e) from the independent-pixel retrieval contains many artifacts (high spikes and low dips in O 3 columns), which are substantially reduced using the SOE method, resulting in a much more realistic (smooth) O 3 map (Fig. 14f). The O 3 differences ( ) between optimized and independent-pixel retrievals (see Fig. 14c, d, and h) illustrate the quantitative improvements, with a small mean O 3 difference (mean within ±0.5 DU) and a sizable reduction in O 3 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 r s f c extrapolations.
In summary, the SOE method performs single-band (B1 or B2) multiple-IFOV (or image-based) retrieval, yielding reflectivity (R) and total O 3 column ( ), with the associated profile determined by the O 3 model (Eq. 12) that retains only the column-dependent M2TCO3 climatological profile X m . The a priori knowledge of the O 3 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.

Total O 3 column
Radiance measurements of EPIC B1 and B2 radiances have O 3 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 O 3 -sensitive band, both bands jointly provide more information that allows the refinement of climatological representation of the O 3 profile. This refinement is performed by adjusting the climatological profile with the most probable modal profile (e 1 , see Eq. 12) so that both B1 and B2 yield the same total O 3 column.
For retrieval from EPIC, the full state vector to be inverted is x = { 0 , γ 1 , , R 1 , R 2 }, where 0 is the total O 3 column, γ 1 the O 3 profile adjustment factor, the total vertical SO 2 column, and R 1 and R 2 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 O 3 vertical column is estimated first assuming there is no SO 2 . The iteration starts with an initial state vector where c is the climatological total column selected from the M2TCO3 climatology based on time and location. R S 1 and R S 2 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(σ 2 B1 = 0.00345 2 , σ 2 B2 = 0.00345 2 ), as estimated from the random errors of the radiance (I M ) measurements (see Sect. 6.2).
There is no correlation among retrieval parameters: total O 3 column ( 0 ), the deviation (ω 1 ) of the O 3 profile from the mean, SO 2 column ( ), and the MLER parameters R, except between R 1 and R 2 . The diagonal elements of the a priori covariance matrix are S a = diag(ε 2 0 = 10 2 DU 2 , ε 2 γ 1 = 2 2 DU 2 , ε 2 = 0.0001 2 DU 2 , ε 2 R1 = 0.001 2 , ε 2 R2 = 0.001 2 ). The off-diagonal elements are equal to zero, {S a (i, j ) = 0, when i = j }, except for the elements associated with R 1 and R 2 , which may be set at S a (4, 5) = S a (5, 4) = 0.98 ε R1 ε R2 = 0.0099 2 , representing a high degree of correlation (0.99) between R 1 and R 2 . This S a essentially limits the adjustments at each iteration: | 0 | ε 0 (10 DU), | ω 1 | ε γ 1 (2 DU), | | ε (10 −4 DU), | R 1 | ε R1 (0.001), and | R 2 | ε R2 (0.001). The strong constraint on SO 2 ensures that its column never deviates far (> 0.01 DU) from its initial value of 0 during the iteration, essentially enforcing an SO 2 -free retrieval. The strong constraints on R 1 and R 2 also ensure that they remain nearly the same as their initial values R S 1 and R S 2 . The constraints on O 3 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 S a .
With the setup of error and a priori covariance matrices S and S a , the initial state vector x 0 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 O 3 column ( ) is obtained by integrating the profile X = X m ( 0 )+γ 1 e 1 ( 0 ). In processing EPIC data, the initial O 3 column c of x 0 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 O 3 , LER, and AI, which are respectively shown in Fig. 15c, e, and f. For comparison, the MERRA-2 assimilated O 3 total columns, interpolated to the time and location of EPIC IFOVs, are included in Fig. 15d, their differences O 3 (EPIC) − O 3 (MERRA-2) in Fig. 15g, and the histogram in Fig. 15h. This comparison reveals excellent agreement between MERRA-2 and EPIC total O 3 , showing nearly identical O 3 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 O 3 processing (see the procedure in Algorithm 1), showing the total O 3 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 O 3 fields, with difference statistics showing slightly worse offsets (µ(EPIC B1) = 0.25 % and µ(EPIC B2) = −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 O 3 spread by σ 2 (EPIC B1) − σ 2 (EPIC) = 0.9 % (or 2.8 DU) and the B2 O 3 spread by a similar amount. These better agreements are consistent over time and location, substantiating the improved retrieval with both O 3sensitive bands over a single one, which is adopted by the TOMS-V8 algorithm. Since the MERRA-2 O 3 field from the assimilation of independent measurements of the Aura OMI and Aura Microwave Limb Sounder (MLS) provides highly realistic spatiotemporal O 3 representation, the smaller spread between the two-band (B1 and B2) EPIC and MERRA-2 total O 3 columns indicates that the inclusion of more O 3sensitive bands enables more accurate retrievals.

Volcanic SO 2
EPIC B1 and B2 radiances respond to both O 3 and SO 2 absorptions, but with very different (see Fig. 13) sensitivities: SO 2 is more than twice as UV-absorbent as O 3 at B1; in contrast, it is significantly less at B2 and about 70 % as absorbent as O 3 . Consequently, the estimate of O 3 absorption signals at these two bands would result in an error due to the presence of SO 2 in the atmosphere: 1 DU of SO 2 would usually yield more than 2 DU O 3 error at B1, but only about 0.7 DU error at B2. This big difference in absorption sen-sitivities facilitates the detection of SO 2 in the atmosphere. Given a radiance signal-to-noise ratio (SNR) of 290 : 1, the theoretical minimum detectable level of SO 2 enhancement is ∼ 0.5 DU in the upper troposphere and above. However, it is difficult to distinguish SO 2 at this minimum level from other changes, such as the O 3 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 SO 2 in the atmosphere. Consequently, low levels of SO 2 elevation cannot be reliably detected in EPIC observations. For significant SO 2 elevations, typically those from volcanic eruptions, B1 O 3 is much higher than B2 O 3 (i.e., 1 > 2 ) from the total O 3 retrieval (described in Sect. 5.2). Adjusting the O 3 profile shape or changing the spectral reflectance of the underlying surface usually cannot eliminate this large O 3 discrepancy between the two bands. Therefore, a high positive value of can be used to flag the presence of SO 2 . Furthermore, a volcanic plume usually occupies a contiguous area with a limited spatial extent. Thus, and 1 enhancements resulting from volcanic SO 2 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 SO 2 plumes, we next describe an algorithmic procedure to flag IFOVs with SO 2 enhancements.
For reliable SO 2 detection, the following procedure is applied to identify the presence of SO 2 in an IFOV. First, IFOVs of likely SO 2 elevations are flagged through spatial analysis of the differential O 3 field (i.e., = 1 − 2 , EPIC B1 and B2 O 3 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 SO 2 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 km 2 that covers the contour). Within this closed 1 contour, IFOVs with likely SO 2 enhancements are flagged when 1 > 2 . This flagging is then extended to the adjacent areas outside the two contours to identify IFOVs with possible SO 2 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 SO 2 elevations.
Once detected, the SO 2 quantification follows the DVCF retrieval with the initial state and a priori covariance setting described next. For the IFOV identified with SO 2 contamination, the initial O 3 values are spatially interpolated from the background 1 field, γ 1 = 0, and initial SO 2 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 S a = diag(ε 2 0 = 10 2 DU 2 , ε 2 γ 1 = 2 2 DU 2 , ε 2 = 0 2 DU 2 ), i.e., the variances associated with O 3 are the same as those for total O 3 retrieval, while the SO 2 column variance is equal to the square of the initial SO 2 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 O 3 retrieval described in the previous section. We list the complete algorithmic procedure in Algorithm 2, titled EPIC total SO 2 retrieval, for SO 2 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 SO 2 product. Figure 16 illustrates the detection and quantification of volcanic SO 2 from EPIC observations. Spatial analyses (i.e., contour mappings) of the intermediate results ( Fig. 16a and b) of the total O 3 processing (see Algorithm 1) provide reliable detection of SO 2 elevations. The SO 2 -flagged IFOVs are then processed with the DVCF algorithm to retrieve total vertical O 3 and SO 2 columns simultaneously, with results shown in Fig. 16c and d, respectively. Comparison of the two O 3 fields in Fig. 16 shows that the initial O 3 elevations (Fig. 16a) due to the presence of SO 2 are nearly entirely removed in the final O 3 field (Fig. 16c), demonstrating that the combined retrieval of O 3 and SO 2 achieves consistent O 3 values inside and outside the plumes. The achieved internal consistency indirectly validates the SO 2 columns. In Fig. 16e, we show maps of DVCF-retrieved SO 2 columns from a series of eight consecutive EPIC observations of the Raikoke plume in Fig. 16e, with the maximum SO 2 value, the total SO 2 mass, and the total area covered by elevated SO 2 displayed in each snapshot. These results illustrate the high-cadence observing capability and high-quality SO 2 measurements of EPIC.

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 O 3 and SO 2 products.

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 (x L ): where the column vectors of length n l (the number of atmospheric layers), ω, and ξ respectively represent the actual vertical profiles for O 3 and SO 2 , while ω L is the climatological O 3 profile equal to X m ( L ) + γ 1 e 1 ( L ) and ξ L the prescribed SO 2 profile specified by a GDF layer with an integrated vertical column equal to L . The profile weighting functions, k ω = −∂ ln I TOA ∂ω | ω=ω L and k ξ = −∂ ln I TOA ∂ξ | ξ =ξ L , are m × n l 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 b L are respectively a set of the exact forward model parameters and those used in the linearization, with the corresponding sensitivity matrix k b = −∂ ln I TOA ∂b | b=b L . 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 O 3 and SO 2 , 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 −x L and putting Eq. (32) into residual y L , Eq. (25) is rewritten as where A ω = Gk ω and A ξ = Gk ξ are the averaging kernels (AKs) for O 3 and SO 2 , respectively. Equation (34)  (c) Vertical O 3 column from EPIC total SO 2 retrieval (see Algorithm 2). (d) Vertical SO 2 column from EPIC total SO 2 retrieval. (e) SO 2 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. and written as after subtracting T − L and T − L from the row equations, respectively. Here T and T are the true O 3 and SO 2 columns, integrated from the corresponding true O 3 (ω) and SO 2 (ξ ) profiles. A and G are the row vectors associated with the retrieved O 3 column from the corresponding matrices A ω and G. Analogously, A and G are the row vec-tors related to the retrieved SO 2 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 O 3 accuracy, including errors in model parameters k b (b − b L ), forward modeling f , spectral measurements m , and the other absorber k ξ (ξ − ξ L ). Similarly, represents the combined total error affecting total SO 2 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 O 3 and SO 2 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 O 3 and SO 2 retrievals contributed from various error sources and simplified physics treatments.

Measurement errors
Errors in EPIC spectral measurements contribute to uncertainties in retrieved O 3 and SO 2 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 O 3 -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 O 3 and SO 2 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 de- Figure 17. Noise levels, i.e., standard deviations (σ ) of O 3 and SO 2 errors ( ), contributed from the random noises in EPIC spectral measurements. The SO 2 noise estimate is for a layer at an altitude 11 km above sea level.
viations, σ ( ) 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 O 3 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 SO 2 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.

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 (m z σ O 3 ,SO 2 ) and (m z ∂σ O 3 ,SO 2 ∂T ) are m × n l matrices with elements {m z j (λ i ) σ O 3 ,SO 2 (λ i , T j ), i = 1. . .m, j = 1. . .n l } and {m z j (λ i ) ∂σ O 3 ,SO 2 (λ i ,T ) ∂T | T j , i = 1. . .m, j = 1. . .n l }, respectively. Here m z (see Eq. 7) is the mean photon path length through a layer at altitude z, σ O 3 ,SO 2 represents errors in O 3 or SO 2 cross-sections, and T represents errors in atmospheric temperature profiles. 5902 X. Huang and K. Yang: O 3 and SO 2 ATBD The BDM O 3 cross-sections (Daumont et al., 1992;Brion et al., 1993;Malicet et al., 1995, commonly abbreviated as BDM) and the BW SO 2 cross-sections (Birk and Wagner, 2018) are used in O 3 and SO 2 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 O 3 cross-sections (Birk and Wagner, 2021) and the SO 2 absorption crosssections of Bogumil et al. (2003) are the alternatives that can replace the baselines for EPIC retrievals. The cross-section errors, σ O 3 in Eq. (39) and σ SO 2 in Eq. (40), are estimated based on the differences between the alternatives and baselines, showing that alternative O 3 and SO 2 crosssections 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 O 3 column biases between 0.5 % and 2 % as well as SO 2 column biases between 1 % and 2 % depending on the effective cross-section differences. The temperature dependence of the BW O 3 cross-sections behaves quite differently from the BDM (Bak et al., 2020), especially at EPIC B2, contributing to the high ends of O 3 biases (> 1.0 %), which occur predominantly at high (viewing or solar) zenith angles when O 3 retrieval becomes more sensitive to EPIC B2 radiance.
Both O 3 and SO 2 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 Liu, 2019) 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 Liu, 2019), 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.

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 (Dave, 1964) with pseudo-spherical approximation (Caudill et al., 1997) to the problem of the transfer of so-lar 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 O 3 and SO 2 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.

Profile errors
As described in Sect. 3, a column-dependent O 3 profile, whose tropospheric integration matches the climatological tropospheric column, is used to specify the vertical distribution of a retrieved total O 3 column. This retrieved profile (ω L ), which represents the likely vertical distribution of the retrieved O 3 vertical column, invariably differs from the actual profile (ω). The O 3 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 O 3 retrievals are shown in Fig. 18a and c, illustrating how A changes with viewing geometry. For low VZAs (< 55 • ), O 3 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, O 3 AKs change with surface reflectance in addition to angular dependence. Under cloud-free conditions, O 3 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), O 3 AKs increase drastically (see Fig. 18c), indicating enhanced sensitivity to tropospheric profile, especially at low zenith angles. Evidently, O 3 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 lowreflectivity surface. In general, errors due to the profile shape are reduced for high-reflectivity surfaces.
Over a short period (e.g., 1 d), O 3 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 O 3 errors ( ) can be written as where the expected values (i.e., the statistical means), E (ω − ω L )(ω − ω L ) T , are O 3 profile covariance matrices S n l , which depend on total columns ( ), season, and latitude. This random component is estimated as a function of VZA using the column-dependent S n l from the M2TCO3 climatology (Yang and Liu, 2019). 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. Retrieval of SO 2 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 SO 2 retrieval error, which can be estimated using the retrieval AK. Figure 18b and d show sample AKs and their variations with VZA for an SO 2 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 SO 2 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.

Errors from Lambertian treatment of natural surfaces
Reflections from surfaces are anisotropic but treated as isotropic. To estimate errors in O 3 and SO 2 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 (Spurr, 2006), for a molecular atmosphere with various O 3 and SO 2 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 O 3 and SO 2 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 O 3 and SO 2 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 O 3 and SO 2 . Test results (see, e.g., Fig. 20) show that errors in total O 3 are mostly within ±1 DU, while SO 2 errors are within ±5 % for SO 2 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 O 3 profile is affected by this approximation. Thus, O 3 errors are proportional to the tropospheric columns but are insensitive to the total column amounts. Since a vast majority of volcanic SO 2 clouds are located below 20 km altitude, SO 2 errors are proportional to the total SO 2 columns. Higher SO 2 clouds are not affected by this treatment.

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 O 3 and SO 2 profiles over Lambertian surfaces of different reflectivities. Then inversion from the simulated radiances with Figure 18. Examples of EPIC total O 3 AKs (a, c) and SO 2 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 O 3 and 30 DU of SO 2 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.
the MLER treatment permits the identification of conditions under which retrieval errors are significant.

Clouds
The error in total O 3 due to the MLER treatment of a lowlying (below 10 km) cloud is mostly within ±2 DU (e.g., Fig. 21a). This O 3 error decreases slightly with a lower cloud altitude (or higher cloud pressure) but is insensitive to the cloud fraction or the total O 3 column. In other words, the MLER treatment does not contribute to large uncertainty in the retrieved total O 3 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 O 3 column: a low (high) bias in OCP results in a positive (negative) error in total O 3 ; 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 Figure 19. The standard deviation (σ ) of O 3 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). the total O 3 uncertainty. Combining O 3 uncertainties due to the OCP error and the MLER treatment yields ±4 DU uncertainty in total O 3 under cloudy conditions. The error in the total SO 2 column due to the MLER cloud treatment is within ±2 % when the SO 2 layer in the troposphere is well above the underlying cloud. This SO 2 error increases with a smaller separation between the SO 2 layer and the cloud, reaching ±15 % when the SO 2 layer is just above the cloud. These characteristics of SO 2 error are illustrated in Fig. 21b. In contrast to the MLER treatment O 3 error, which is insensitive to the total column, this SO 2 error is proportional to the total SO 2 column. When the SO 2 layer is below or within the cloud, the uncertainty of SO 2 quantification increases drastically. Depending on the relative distributions of SO 2 and the cloud particles, the retrieved SO 2 based on the MLER treatment can be a fraction of or a few times the actual column.

Aerosols
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 nonabsorbing and weakly absorbing aerosols in the lower troposphere (< 7 km) results in small (< ±2 DU) errors in total O 3 retrievals, provided that the proper OCP for the elevated cloud surface is used. This error range is nearly independent of the total O 3 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 r g f c extrapolation scheme (see Sect. 2.3) result in a positive bias in the retrieved total O 3 columns (see, e.g., Fig. 22). This O 3 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 O 3 profiles reveals a positive correlation between column O 3 error and the UV AI. Quantitatively, this relationship can be written as O 3 = (1.5 ± 1) × AI DU for AI values greater than 0.5 and less than 8. This relationship provides a rough estimate of O 3 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 O 3 bias of about 3 DU for IFOVs contaminated with UVabsorbing 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.
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 SO 2 ) into the atmosphere. Since ash particles strongly absorb UV, the MLER treatment of volcanic plumes leads to huge uncertainties in the retrieved SO 2 columns, which are often greatly overestimated or underestimated depending on the relative distributions between SO 2 and ash particles. An explicit treatment of volcanic ash is needed for accurate retrieval of SO 2 when ash particles are co-located with or slightly separated from the gas.

Error summary
Uncertainty estimates (Sect. 6.2) have detailed various contributions to systematic and random errors in retrieved O 3 and SO 2 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 O 3 vertical profiles, are combined to estimate total random errors. The random O 3 error (characterized by its standard deviation) is 1.5 % at low VZAs (< 45 • ), increasing to 2.0 % at 75 •  . Errors in retrieved O 3 and SO 2 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) O 3 errors in DU for correct (p c = 545 hPa) and biased (p c = 545 ± 100 hPa) cloud OCPs. (b) SO 2 errors in percent for SO 2 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.
for IFOVs without clouds and aerosols. Similarly for IFOVs with clouds or non-absorbing aerosols, the random error in the total O 3 column is 1.8 % for low VZAs (< 45 • ) and 2.5 % at 75 • . Random SO 2 error has two terms: one is independent of the SO 2 column (see Sect. 6.2.1), but the other is proportional to this column (see Sect. 6.2.2). Hence, they are not combined.
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 O 3 column to be ±2 % for VZA ≤ 70 • and ±3.5 % for VZA ≤ 85 • ; the bias in the total SO 2 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, O 3 has a positive bias for IFOVs with UV-absorbing aerosols, increasing with a higher AI (see Sect. 6.2.6). Under-estimation (overestimation) of an SO 2 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 O 3 and comparison of SO 2

O 3 validation
We validate the DVCF O 3 retrievals from EPIC using ground-based Brewer spectrophotometer measurements and the assimilated O 3 product from MERRA-2, the Modern-Era Retrospective Analysis for Research and Applications Version 2 (Gelaro et al., 2017).
We compare EPIC total O 3 columns with the Brewer O 3 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 O 3 columns that are coincident (within ±15 min) with EPIC observations at these stations. For intercomparison, Brewer O 3 data are interpo- lated to the times when EPIC observes these locations. Coincident O 3 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 O 3 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 O 3 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 O 3 variability accurately. EPIC O 3 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 O 3 .
From October 2004, MERRA-2 O 3 field is assimilated from Aura MLS and OMI and provides highly realistic global distributions of O 3 in the stratosphere and upper troposphere while inheriting the uncertainty characteristics of its sources (Stajner et al., 2008;Davis et al., 2017). We compare the MERRA-2 synoptic O 3 field with the EPIC hemispheric view for the same observation time to access EPIC's capability of capturing the realistic O 3 distribution. For instance, strikingly similar O 3 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 O 3 , we estimate the accuracy of EPIC total O 3 to be −0.76 ± 2.19 %.

SO 2 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 SO 2 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 polarorbiting 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 O 3 and SO 2 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 SO 2 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 SO 2 quantification because the SO 2 measurement sensitivity varies strongly with its altitude. Thus, DVCF height-resolved SO 2 retrievals from hyperspectral instruments, such as OMI and OMPS-NM, provide the most accurate quantification of SO 2 vertical columns. To validate DVCF SO 2 retrievals from EPIC, we compare the SO 2 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 SO 2 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 SO 2 into sulfate, and the changing measurement sensitivity with viewing illumination conditions since low SO 2 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 SO 2 plume, thus usually providing more accurate estimates of the lower bound of SO 2 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 SO 2 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.

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 O 3 , SO 2 , and NO 2 , from Aura OMI and Suomi NPP as well as NOAA-20 OMPS. Algorithm advances, including the improved O 3 profile representation and the regulated direct fitting inversion technique, improve the accuracy of O 3 and SO 2 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 O 3 ) from that (e.g., tropospheric O 3 ) with higher spatial variations.
A thorough error analysis is provided to quantify O 3 and SO 2 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 O 3 and SO 2 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 O 3 and SO 2 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 O 3 columns are validated against coincident ground-based Brewer measurements and compared with coincident O 3 data from MERRA-2 assimilation. The findings show that EPIC total O 3 is highly accurate, realistically capturing the short-term O 3 variability while maintaining long-term consistency over the entire record. The EPIC SO 2 loadings of volcanic plumes are evaluated against those from hyperspectral measurements of the same eruptions, showing that EPIC provides accurate SO 2 quantifications from large volcanic eruptions. EPIC's high-cadence observations allow better identification of the peak loading of volcanic SO 2 plumes compared to polar-orbiting instruments.
The BW O 3 and SO 2 cross-sections Wagner, 2021, 2018)  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 algo-