A new discrete wavelength backscattered ultraviolet algorithm for consistent volcanic SO2 retrievals from multiple satellite missions

Abstract. This paper describes a new discrete wavelength algorithm
developed for retrieving volcanic sulfur dioxide (SO2) vertical column
density (VCD) from UV observing satellites. The Multi-Satellite SO2
algorithm (MS_SO2) simultaneously retrieves column densities
of sulfur dioxide, ozone, and Lambertian effective reflectivity (LER) and its
spectral dependence. It is used operationally to process measurements from
the heritage Total Ozone Mapping Spectrometer (TOMS) onboard NASA's
Nimbus-7 satellite (N7/TOMS: 1978–1993) and from the current Earth
Polychromatic Imaging Camera (EPIC) onboard Deep Space Climate Observatory
(DSCOVR: 2015–ongoing) from the Earth–Sun Lagrange (L1) orbit. Results from
MS_SO2 algorithm for several volcanic cases were assessed
using the more sensitive principal component analysis (PCA) algorithm. The
PCA is an operational algorithm used by NASA to retrieve SO2 from
hyperspectral UV spectrometers, such as the Ozone Monitoring Instrument (OMI) onboard NASA's Earth Observing System Aura satellite and Ozone Mapping and
Profiling Suite (OMPS) onboard NASA–NOAA Suomi National Polar Partnership
(SNPP) satellite. For this comparative study, the PCA algorithm was
modified to use the discrete wavelengths of the Nimbus-7/TOMS instrument,
described in Sect. S1 of the Supplement. Our results demonstrate good
agreement between the two retrievals for the largest volcanic eruptions of
the satellite era, such as the 1991 Pinatubo eruption. To estimate SO2
retrieval systematic uncertainties, we use radiative transfer simulations
explicitly accounting for volcanic sulfate and ash aerosols. Our results
suggest that the discrete-wavelength MS_SO2 algorithm,
although less sensitive than hyperspectral PCA algorithm, can be adapted to
retrieve volcanic SO2 VCDs from contemporary hyperspectral UV
instruments, such as OMI and OMPS, to create consistent, multi-satellite,
long-term volcanic SO2 climate data records.



Introduction
Volcanic eruptions are an important natural driver of global climate change, but, unlike other natural climate forcing (e.g., changes in Earth's orbit, solar irradiance), the magnitude of volcanic forcing is highly variable and largely unpredictable, and the effects are typically more transient. Of most interest are the episodic, large injections of volcanic sulfur dioxide (SO 2 ) into the Earth's stratosphere by major explosive volcanic eruptions, the most recent example being the eruption of Pinatubo (Philippines) in June 1991 (e.g., Bluth et al., 1992;Guo et al., 2004). Stratospheric loading of volcanic SO 2 by major eruptions leads to the formation of sulfuric acid (or sulfate) aerosols that scatter incoming solar shortwave radiation and absorb outgoing thermal radiation over timescales of months to years, cooling the troposphere and warming the stratosphere (e.g., Robock, 2000). Primary volcanic emissions of aerosols such as volcanic ash can also have atmospheric and climatic impacts, but these are typically more short-lived. Volcanic eruptions can also release reactive halogen species into the atmosphere, such as chloride and bromide (Mankin and Coffey, 1984;Bobrowski et al., 2003;Kern et al., 2009). Halogens can impact the total column ozone amount and profile shape if injected into the lower stratosphere (Solomon et al., 1998;Klobas et al., 2017), but sulfate aerosols are also required to catalyze the heterogeneous chemical reactions that can efficiently deplete ozone. Hence, to understand the impacts of volcanic eruptions on climate, and in order to predict possible outcomes in the event of a major eruption, long-term satellite measurements of volcanic SO 2 emissions are essential.
The satellite record of volcanic SO 2 emissions by major volcanic eruptions extends back to 1978 and has been derived from instruments operating in both the ultraviolet (UV) and infrared (IR) spectral bands ( Fig. 1; e.g., Carn et al., 2003Carn et al., , 2016Carn, 2019;Prata et al., 2003). Measurements in the UV band have a longer heritage, since the first satellite detection of volcanic SO 2 was achieved by the UV Total Ozone Mapping Spectrometer (TOMS) in 1982 following the eruption of El Chichón (Mexico; Krueger, 1983;Krueger et al., 2008), and interference from volcanic SO 2 must be accounted for in order to produce accurate, long-term UV measurements of ozone. UV measurements have greater sensitivity to the total atmospheric SO 2 column than IR retrievals, and hence the former have been the mainstay of volcanic SO 2 monitoring during the satellite era to date. The volcanic SO 2 climatology from 1978 to the present (Fig. 1, Carn, 2019) reveals highly variable inter-annual volcanic SO 2 forcing dominated by two major eruptions (El Chichón in 1982 andPinatubo in 1991), with the post-2000 period dominated by smaller eruptions. Although none of these smaller eruptions have, individually, produced measurable climate effects, collectively they have garnered significant interest as they may play an important role in sustaining the persistent, background stratospheric aerosol layer, which is an important factor in global climate forcing (e.g., Solomon et al., 2011;Vernier et al., 2011;Ridley et al., 2014).
One of the key challenges in assembling a long-term, consistent, satellite-based volcanic SO 2 emissions climatology (e.g., Fig. 1) is merging measurements from sensors with different spectral coverage and resolution. This complicates any analysis of "trends" in volcanic SO 2 loading (e.g., in the post-2000 period of moderate volcanic eruptions; Fig. 1) or comparisons of eruptions of similar magnitude in different satellite instrumental eras. A step change in SO 2 sensitivity occurred when the multispectral, six-channel TOMS instruments were superseded by hyperspectral UV sensors, such as the Global Ozone Monitoring Experiment (GOME, 1995(GOME, -2003Khokhar et al., 2005), the Scanning Imaging Absorption Spectrometer for Atmospheric Chartography (SCIA-MACHY, 2002Lee et al., 2008), the Ozone Monitoring Instrument (OMI, 2004-ongoing;Krotkov et al., Figure 1. Multi-decadal record of SO 2 emissions by volcanic eruptions observed by NASA's fleet of satellites observing TOA UV radiances. Eruptions (star symbols) are color coded by estimated plume altitude and derived from a variety of sources, including Smithsonian Institution Global Volcanism Program volcanic activity reports, volcanic ash advisories and satellite data. The annual total explosive volcanic SO 2 production (omitting SO 2 discharge from effusive eruptions) is shown in black. The orange lines above the plot indicate the operational lifetimes of NASA UV satellite instruments: Nimbus-7 (N7), Meteor-3 (M3), ADEOS (AD), Earth Probe (EP) TOMS, OMI (currently operational), and SNPP/OMPS (currently operational), along with the ESA/EU Copernicus S5P/TROPOMI (currently operational). Data shown in this plot are available from the NASA Goddard Earth Sciences (GES) Data and Information Services Center (DISC) as a level 4 MEaSUREs (Making Earth System Data Records for Use in Research Environments) data product (Carn, 2019). 2006), the Ozone Mapping and Profiler Suite (OMPS, 2012ongoing;Carn et al., 2015), and EU/ESA Copernicus Sentinel 5 precursor (S5P) (Veefkind et al., 2012). This is manifested in Fig. 1 as an increased number of detected volcanic eruptions with low SO 2 loading (< 10 kt) after 2004 (note that GOME and SCIAMACHY measurements are not shown in Fig. 1), whereas rates of global volcanic activity have not changed significantly. UV SO 2 retrieval algorithms have also evolved substantially since the 1980s in response to advances from multispectral to hyperspectral sensors, improvements in ozone retrievals, and efforts to account for volcanic ash and aerosol interference (e.g., Krueger et al., 1995Krueger et al., , 2000Krotkov et al., 1997Krotkov et al., , 2006Yang et al., 2007Yang et al., , 2010Li et al., 2013Theys et al., 2015). However, to date there has been no attempt to develop a single algorithm that could be used to generate a long-term, consistent SO 2 climatology across multiple UV satellite missions. In this paper we describe a new multi-satellite SO 2 algorithm (MS_SO2) that is applicable to both multispectral (e.g., TOMS) and hyperspectral (e.g., OMI) UV measurements. As a first step in the generation of a multi-satellite volcanic SO 2 record, we apply the MS_SO2 algorithm to the Nimbus-7 TOMS (N7/TOMS) measurements (1978)(1979)(1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993) and present a reanalysis of some of the most significant eruptions of the N7/TOMS mission.

Heritage satellite ozone and SO algorithms
Ozone and SO 2 are the two main absorbers in the near UV spectral region between 300 and 340 nm. The relative contributions of each gas to the satellite backscattered ultraviolet (BUV) measurements at the three absorbing TOMS channels (317, 331, 340 nm) used in the retrieval, depend on the spectral structure of the absorption cross sections, which are measured as functions of wavelength and temperature (Bogumil et al., 2003;Daumont et al., 1992). Figure 2 shows the O 3 and SO 2 cross sections and the SO 2 /O 3 cross section ratio as a function of wavelength for a spectral UV region spanning the three absorbing channels of TOMS. At the instrument's spectral resolution (∼ 1 nm full width at half maximum, FWHM), the SO 2 molecule is 2.5 times more absorbing than O 3 at 317 nm, while O 3 is 6 times more absorbing at 331 nm. These differences allow for simultaneous multispectral retrievals of O 3 and SO 2 . Dave and Mateer (1967) first proposed a technique to estimate total ozone column from nadir backscatter UV measurements taken in the Huggins ozone absorption band (310-340 nm), assuming no SO 2 is present. Their algorithm was inspired by the pioneering Dobson spectrophotometer, which measures attenuation of solar irradiance by UV wavelength pairs from which total ozone is derived using the Beer-Lambert law. However, unlike the direct sun technique, radiative transfer calculations show that the top-of-the-atmosphere (TOA) BUV radiances (I ) do not follow the Beer-Lambert law. In general, log(I ) varies nonlinearly with ozone column amount ( ), and this relationship is sensitive to the shape of the ozone profile (defined as the ozone density profile normalized to total ozone). To account for this effect, Dave and Mateer (1967) proposed constructing a set of lookup tables (LUTs) based on standard ozone profiles with different total ozone amounts using ozonesonde and Dobson Umkehr data. Since the shapes of the profiles also vary with latitude, they proposed using three sets of profiles for low latitudes, midlatitudes and high latitudes. These profiles are then used to estimate I , which varies with wavelength (λ), observational geometry, surface pressure and surface reflectivity (R). Following the Dobson convention, log(I ) is converted to N value which is defined in Eq. (1) as

Heritage BUV ozone algorithms
(1) F is the extraterrestrial solar irradiance. By linearly interpolating N between total ozone nodes, one forms the Ncurves that are a single valued function of representative of a given latitude band and observational geometry. This approach allows to be estimated by matching the measured N value to the interpolated N values. Over the years several modifications have been introduced to this basic concept. Mateer et al. (1971) proposed a Lambert-equivalent reflectivity (LER) concept to estimate the combined contribution of the surface, clouds and aerosols to BUV radiance. In this concept, the scene at the bottom of the atmosphere is assumed to be a Lambertian reflector, whose reflectivity (R s ) is derived from the measurements at 380 nm where the ozone and SO 2 absorption is negligible. The effective pressure of this reflecting surface is assumed to vary with R s , from a surface pressure at R s < 0.2 to a cloud pressure 0.4 atm at R s > 0.6, linearly interpolated at intermediate R s . The algorithm assumed that R s , thus derived, did not vary with wavelength. Although in the earlier versions of this algorithm wavelength pairs (313/331, 318/340) were used to derive , R s was later derived at 331 nm to minimize errors due to the spectral dependence of R s . This made pairing unnecessary (McPeters et al., 1996).
By explicitly modeling the effect of aerosols using a radiative transfer code, Dave (1978) showed that R s did not vary significantly with wavelength for non-absorbing aerosols, hence they produced no ozone error. However, for aerosols that might have strong absorption in the UV, that study predicted that R s would decrease at shorter wavelengths, producing an overestimation of ozone. Since aerosol properties in the UV were not known at that time, no correction for aerosol absorption was applied until the mid-1990s when the effect predicted by Dave (1978) was detected in the Nimbus-7 TOMS data launched in October 1978.
Since the TOMS instrument had three reflectivity channels (331, 340, 380 nm), it was possible to compare the reflectivities derived from them. This comparison showed that R s in-creased significantly with wavelength for moderately thick clouds, causing a significant underestimation of (up to 3 %). A modified LER (MLER) concept assuming two Lambertian surfaces, one at the surface and the other at the cloud top was applied to minimize this error (Ahmad et al., 2004).
The most recent version of the TOMS ozone algorithm reverts back to the LER model, but it assumes that clouds are at the surface, which reduces the R s wavelength dependence (Ahmad et al., 2004). This simple LER (SLER) model is used in our SO 2 algorithm. However, since there are many other reasons for such a dependence, including ocean color, non-Lambertian surfaces, such as ocean glint and fogbow and, most importantly, the absorbing aerosol effect predicted by Dave (1978), R s is assumed to vary linearly with λ; its slope is derived using 340 and 380 nm radiances. This simple omnibus approach works well for most cases, except when the UV absorbing aerosols (smoke, dust and volcanic ash) are very thick. Such data are flagged in the TOMS ozone algorithm. The new MS_SO2 algorithm is an extension of this algorithm into two dimensions (Sect. 3). Krueger (1983) was the first to suggest that TOMS could be used to retrieve sulfur dioxide from explosive volcanic events. He correctly interpreted the large positive ozone anomaly observed following the explosive eruption of El Chichón in 1982 as being due to the SO 2 released into the atmosphere during the event. To estimate the SO 2 inside the plume region, he separated the SO 2 and O 3 signals by computing a residual reflectance, estimated as the difference between the interpolated unperturbed background reflectances outside the plume and the reflectance anomaly inside the plume. This early technique for retrieving SO 2 from TOMS ozone estimates became known as the residual method. The residual method, however, failed when the background could not be clearly separated from the ozone anomaly. Krueger subsequently developed the first BUV algorithm that separated the O 3 and SO 2 radiance contributions, based on an earlier methodology developed by Kerr et al. (1980) to retrieve the SO 2 column from the ground with a Brewer spectrophotometer. This method assumed that the BUV radiation was attenuated by the two absorbing species (O 3 , SO 2 ), leading to an equation describing BUV radiance, I , for a given wavelength, λ, corresponding to the TOMS field of view (FoV):

Heritage TOMS SO 2 algorithms
In Eq.
(2), F is the incoming solar flux, S g is the geometrical optical path (air mass factor, AMF), and τ O 3 and τ SO 2 are the vertical optical thicknesses for O 3 and SO 2 , while the coefficients a and b depend on the satellite viewing geometry, cloud or surface reflectance, and volcanic ash and sulfate aerosols (Krueger et al., 1995;Krotkov et al., 1997). Equation (2) can be expressed in matrix form, which is then inverted to obtain estimates for the SO 2 and O 3 vertical col-umn densities and the dimensionless parameters a and b. This algorithm is generally referred to as the Krueger-Kerr algorithm (Krueger et al., 1995). Krotkov et al. (1997) developed radiative transfer path correction, which explicitly accounted for the R s , ozone and SO 2 vertical profiles, replacing the geometrical AMF in Eq.
(2). The modified algorithm with empirical background correction has been used offline on a case-by-case basis for the past 2 decades to retrieve SO 2 mass tonnage from medium to large explosive eruptions using TOMS BUV measurements Carn et al., 2003).

New MS_SO2 algorithm
The new discrete wavelength SO 2 algorithm (MS_SO2) builds on the heritage of the TOMS total ozone algorithm (Sect. 2.1) but adds sulfur dioxide (SO 2 ) as a second absorber. The BUV radiance is simulated with the TOMRAD forward vector radiative transfer (RT) model (Dave, 1964) for a known viewing geometry by assuming a vertically inhomogeneous, pseudo-spherical Rayleigh scattering atmosphere with standard ozone profiles (Klenk et al., 1983;Bhartia, 2002) and a priori SO 2 vertical profiles (Krueger et al., 1995). The underlying reflecting surfaces (land/ocean, clouds and aerosols) are approximated with the simple LER reflecting surface at terrain height pressure (Sect. 2.1). TOM-RAD accounts for all orders of polarized Rayleigh scattering and for the gaseous absorption (e.g., O 3 and SO 2 ), using a priori vertical profiles of the gas concentrations and laboratory-measured temperature-dependent gaseous cross sections (Dave and Mateer, 1967;Bogumil et al., 2003;Daumont et al., 1992). Improvements to the TOMRAD model include corrections for molecular anisotropy (Ahmad and Bhartia, 1995), rotational Raman scattering (Joiner et al., 1995) and pseudo-spherical corrections to account for changes to the solar and viewing zenith angles due to the sphericity of the Earth. Performing online radiative transfer calculations for every satellite field-of-view (FoV) can greatly increase the time required to process full orbits of data. To improve the computational efficiency of the operational algorithm, N7TOMSspecific lookup tables (N7TOMS-LUT) were produced offline using the inputs listed in Table 1 and convolved with the triangular band pass at each of the six Nimbus-7 TOMS wavelengths (FWHM ∼ 1 nm).
The TOMRAD algorithm was configured to account for two absorbing trace gases: O 3 and SO 2 . The LUTs include 21 total ozone nodes and 12 total SO 2 nodes for each of the three assumed SO 2 heights. For ozone, the total column amounts and profile shapes vary between latitude bands (see Table 1). For sulfur dioxide, we assumed a Gaussian vertical profile shape, which is determined by two parameters: a center of mass altitude (CMA) and a geometrical standard deviation. The CMA represents the altitude of the peak SO 2 where is the retrieved total column sulfur dioxide, is the total column ozone, dR/dλ characterizes the R s spectral dependence between 340 and 380 nm, and R s is the LER at 380 nm. The retrieval of sulfur dioxide is carried out in one or two steps described in the next sections, referred to as step 1 and step 2.

Step 1 retrieval
Our step 1 inversion starts with an initial state vector x 0 , consisting of first guesses for 0 , 0 and dR s0 /dλ shown in Table 2. The final state vector, x, is determined iteratively by inverting the Jacobian matrix K at each iteration step: where dx represents the relative changes in the state vector from the previous iteration and dN represents the residual vector equal to N m −N c , computed as the difference between the measured N values, N m , and the calculated N values, N c , at the four TOMS channels at 317, 331, 340 and 380 nm. K represents a 4 × 4 Jacobian matrix computed from the LUTs. These matrix elements are defined as follows: where N c,i is the forward model calculated N value at wavelength i. The reflectivity R s is computed analytically using the measured BUV radiance at 380 nm (see Supplement, Eq. S4). Note that since the O 3 and SO 2 cross sections are negligible at 380 nm, the R s and ∂N 380 /∂R s do not change with the iterations (i.e., dR s = 0).
Equation (4) is solved iteratively by zeroing the residuals, dN = N m -N c , and recomputing the N c and the Jacobians at each iteration step for the four used channels. The state vector is then adjusted after each iteration, x k = x k−1 + dx k , k = 1, 2,. . . until it converges on a solution as described below: Since O 3 and SO 2 exhibit small absorption at 340 nm, a nonzero R spectral slope (i.e., dR/dλ = 0) accounts for the radiative effects of aerosols and surface reflectance (e.g., sun glint). As indicated in Table 2, the algorithm initially assumes zero R-λ dependence (i.e., dR/dλ = 0); however, absorbing aerosols (smoke, dust and volcanic ash) cause dR/dλ = 0.
The algorithm uses retrieved spectral slope dR/dλ in Eq. (7) below to update the calculated LERs after each iteration: where λ j = 312, 317, 331, 340 nm and λ R = 380 nm. When SO 2 or aerosol loading is high nonlinear R-λ dependence can cause systematic errors in the retrieval state vector. For this reason, we do not use the shortest 312 nm channel in the retrievals (Eqs. 5-6), but the final residual dN 312 is used as a diagnostic of the nonlinearity. A step 2 empirical procedure, described in the next section, was developed to correct for the retrieval bias resulting from these errors.

Step 2 retrieval
The MS_SO2 forward model accounts for O 3 and SO 2 absorption and linear spectral changes in R s due to the presences of aerosols. The algorithm, however, does not explicitly characterize the absorption and scattering effects of volcanic ash (absorbing) and sulfate (non-absorbing) aerosols. The retrieval errors in and caused by volcanic ash during the first few days after an explosive eruption can be significant in the case of major volcanic eruptions like Pinatubo and El Chichón (Krueger et al., 1995;Krotkov et al., 1997). A step 2 procedure was developed primarily to handle explosive eruptions (Volcanic Explosivity Index, VEI > 3 ), in which large anomalies are identified to occur in conjunction with high ash concentrations. In step 2, a corrected total ozone cor inside the SO 2 cloud is interpolated using the retrieved outside the plume along the orbit for each cross-track position. Even if ozone-destroying chemicals are present, such effects can still be considered negligible over the relatively short time periods that SO 2 concentrations are high enough to affect TOMS observations.
In deciding whether to apply step 2, the algorithm considers the retrieved , and aerosol index (AI) in Step 1. The AI is estimated from the dR/dλ and the calculated Jacobian dN/dR at 340 nm: Positive AI (dR/dλ > 0) identifies spatial regions affected by absorbing aerosols (dust, smoke and ash). The step 2 selection criteria first select FoVs where either SO 2 > 15 DU (inside the plume) or AI > 6. The additional AI criterion allows for the selection of FoVs around the edges of the cloud, where the SO 2 can be less than 15 DU due to high aerosol concentrations. In this case, it is assumed that the step 1 SO 2 may have been underestimated due to the ozone error caused by high aerosol concentrations (in these cases, the SO 2 retrieved in step 2 may still not exceed 15 DU and therefore would be excluded from the plume in subsequent mass calculations). We describe the methodology for interpolating cor in Eqs. (S5)-(S7). A second retrieval of SO 2 and dR/dλ is then performed by inversion using the measured 317 and 340 nm radiances while treating the ozone cor as a constant. This constraint on the ozone bounds the SO 2 Jacobians computed from the forward model LUTs. The operational MS_SO2 product files include a step 2 algorithm flag (not applied = 0, applied = 1). To illustrate the effects of the step 2 procedure, we consider the 1982 explosive eruption of El Chichón, which emitted ∼ 7 Tg SO 2 (Krueger et al., 2008) the second largest observed in the satellite era (Fig. 1). Figure 3 shows the retrieved AI map during TOMS overpass of the volcano on 4 April 1982, while it was still erupting. High AI values exceeding a value of 10 correspond to biased high step 1 ozone values (Fig. 4a) and underestimated values (Fig. 4c). Figure 4b shows the step 2 corrected cor , making it consistent with the field outside of the volcanic cloud. Figure 4d shows the step 2 , which is much higher than step 1. As can be seen in this particular example, the step 2 correction can significantly increase the SO 2 mass. In this case, the SO 2 mass increased from 2475 (step 1) to 3637 kt (step 2). Peak values increased from 396 to 549 DU in the aerosol affected region. The biases, d and d , for this case are shown in Fig. S1 in the Supplement.
Step 2 was developed primarily to handle extreme eruptions (VEI > 3), such as El Chichón and Pinatubo, where large anomalies sometimes occur in conjunction with high ash concentrations. In practice, step 2 corrections tend to be small (or nonexistent) for most of the eruptions detected observed during the observation period covered by TOMS.
The corrected step 2 values inside the volcanic cloud shown in Fig. 4b appear to be fairly consistent with the regional unperturbed ozone field, but it should be noted that a few remaining high values in the boundary of the plume still exist, which were not selected for step 2 (Fig. 4b). These pixels were not corrected because the threshold criteria were not met, thus may be underestimated. However, their contribution to the total SO 2 cloud mass is insignificant.
Step 2 follows a methodology similar to the original residual method developed by Krueger (1983), which separated the O 3 and SO 2 contributions by subtracting the measured BUV reflectance in the unperturbed region from the BUV radiance anomaly associated with the SO 2 cloud. The MS_SO2 algorithm corrects the overestimated step 1 ozone inside the plume by correcting the positive ozone bias. Our step 2 procedure is typically only applied when the ash and/or SO 2 loading causes the reflectivity dependence to become nonlinear, as the forward model does not explicitly account for volcanic aerosol absorption. This scenario typically lasts for about 1-3 d following a major explosive eruption, during which total retrieved SO 2 mass is likely to be underestimated and in some cases could even increase with time due to ash and ice fallout and plume dispersion. For such extreme cases we recommend estimating SO 2 to sulfate conversion e-folding lifetime using weeks of measurements of the total SO 2 cloud daily mass and extrapolating it back in time to estimate total SO 2 mass emitted on eruption day. This "day one" time extrapolated SO 2 mass is typically larger than retrieved on days immediately following the eruption .

Soft calibration: N value bias correction
We assume that the background sulfur dioxide is below TOMS detection limit in regions of the atmosphere far away from SO 2 sources (e.g., volcanic, anthropogenic). Random errors associated with the retrieval process, however, are normally distributed around zero. We expect that the true volcanic SO 2 , true and the mean of the distribution, clean , to equal zero such that true = clean = 0.
We examined a sample of 90 TOMS orbits in clean regions of the central Pacific Ocean and found a positive bias of about 3 DU (i.e., clean ∼ 3 DU, Fig. 5). A soft calibration procedure was developed for correcting this bias by applying a small constant N 340 value adjustment to the measured 340 nm BUV radiances. The details of this procedure are described in Sect. S3.3 in the Supplement. Figure 5 shows probability density functions (PDFs) of the step 1 SO 2 before (dashed) and after (solid) applying the correction for 11 November 1981. The mean bias is reduced to < 1 DU after applying the correction.

Random errors and SO 2 detection limit
The random errors in the MS_SO2 retrieval were estimated from the standard deviation in the SO 2 from a large data sample that included 90 central Pacific orbits, spanning a 10-year period between 1980 and 1990. Data were restricted to values between −20 and 20 DU (Fig. 6a). Standard devia-  tions were then computed as a function of the TOMS swath position, as shown in Fig. 6b. Figure 6b can be used to characterize the SO 2 detection limits for TOMS. In this section, we compare the TOMS error distribution with the UV Ozone Mapping Profile Suite Nadir Mapper (OMPS-NM), a hyperspectral UV instrument onboard the Suomi National Polarorbiting Partnership (NPP) and NOAA 20 satellites. For this comparison, we selected 1 month of NPP/OMPS spectral data (central Pacific) and applied the MS_SO2 algorithm using the same four wavelength bands of TOMS (Table 2), which were first convolved with the TOMS bandpass function. Figure 6b shows that TOMS retrieval noise depends on the swath position, varying from ∼ 6 DU at nadir to ∼ 4 DU at higher viewing angles, while OMPS is 2-3 times smaller (∼ 2 DU) and is relatively independent of the cross-track position. Using the MS_SO2 algorithm, we subsequently estimate the SO 2 detection limit for TOMS and OMPS-NM to be about 15 DU and 6 DU (∼ 99 % confidence level), respectively. We note that when applying the Principal Component Algorithm (PCA) (Li et al., 2013) to all the 100-200 wavelengths available from OMPS-NM hyperspectral measurements, the noise is reduced by an order of magnitude to ∼ 0.2-0.5 DU, allowing detection of large anthropogenic point sources (emissions more than ∼ 80 kt yr −1 ) .

Systematic errors in volcanic SO 2 plumes
In this section, we evaluate systematic errors of the MS_SO2 retrievals of volcanic SO 2 . The two most significant errors are caused by volcanic aerosols (ash and sulfate) and incorrect assumptions regarding the SO 2 profile, namely the plume height. The radiance tables used by the algorithm account for ozone and SO 2 absorption but do not account for the absorption and scattering by aerosols. The ash errors can be significant during the first couple days after the initial eruption phase Guo et al., 2004). The precomputed radiance tables used by MS_SO2 assume an SO 2 column amount and an a priori CMA and standard deviation (Sect. 3). An incorrect CMA assumption can cause significant SO 2 errors that vary with viewing geometry, ozone and SO 2 column amounts. We characterize these error sources by applying the MS_SO2 algorithm to synthetic radiances.

Uncertainties due to SO 2 plume height
To understand retrieval errors in MS_SO2 algorithm due to assumed a priori SO 2 profiles, we conducted sensitivity tests using the VLIDORT radiative transfer code for the typical observational conditions in the tropics, midlatitudes and high latitudes. Figure 7 shows column SO 2 Jacobians ∂N/∂ at 317 nm for different SO 2 amounts, , nadir angles and scene reflectance as function of the assumed SO 2 height (center of mass altitude, CMA). The Jacobians generally increase with the CMA, meaning that satellite BUV measurements are more sensitive to SO 2 at higher altitudes. This means that the MS_SO2 algorithm will overestimate (underestimate) the SO 2 column amount if the CMA of the a priori profile is lower (higher) than that of the actual SO 2 profile. On the other hand, the sensitivity of SO 2 Jacobians with respect to CMA is affected by several other factors, particularly SO 2 column amounts, geometry (solar zenith angle and viewing zenith angle), the reflectivity of the underlying surface (R s ) and the CMA itself. In general, the sensitivity of SO 2 Jacobians to CMA is greater for SO 2 plumes with large SO 2 loading (e.g., 300 DU vs. 50 DU) at relatively low altitudes (e.g., CMA of 13 km vs. 18 km), lower reflectivities (e.g., R s of 0.05 vs. 0.50) or at high viewing angles relative to the nadir (e.g., VZA of 60 • vs. 0 • ). For calculations assuming typical midlatitude and high-latitude conditions, we found similar sensitivities of SO 2 Jacobians to CMA. From these calculations, we can estimate the errors in the SO 2 Jacobians at 317 nm, assuming that the standard a priori profiles used in MS_SO2 retrievals (CMA: 13 and 18 km) have a ±2 km error in CMA. The results for the tropics, midlatitudes and high latitudes are summarized in Tables S1, S2 and S3, respectively, in the Supplement. As shown in the tables, for SO 2 plumes from relatively moderate eruptions (∼ 50 DU), the relative errors in SO 2 Jacobians due to the error in the CMA are mostly within ±10 %. But for plumes with large SO 2 loading (∼ 200-300 DU) from explosive eruptions such as Pinatubo, the relative error in SO 2 Jacobians may reach as high as 30 % for pixels near the edge of the swaths that have low reflectivity. Additionally, for pixels with the same reflectivity and VZA, the relative errors due to SO 2 height are greater for midlatitude and high-latitude eruptions than for tropical eruptions.
To quantify the retrieval errors due to inaccuracies in the a priori profiles, we used the top-of-the-atmosphere synthetic radiance data generated by VLIDORT as input to the MS_SO2 algorithm. The retrieved SO 2 and O 3 column amounts were compared with assumed in VLIDORT cal-culations (Tables S4-S7). As shown in the tables, for SO 2 plumes with a modest loading (∼ 50 DU), the relative errors in SO 2 column amounts, due to a 2 km error in the a priori profile are typically 10 % or less, whereas the relative errors in O 3 are within 1 %. For plumes with large SO 2 loadings (200-300 DU), the errors in SO 2 amounts due to a 2 km bias in the a priori profile are typically 5 %-15 % but can reach as high as 30 %-40 % for high-latitude plumes with large SZA and VZA. For extreme conditions at high latitudes (Table S5, 13 km a priori profile vs. 15 km actual profile, SO 2 = 300 DU), the MS_SO2 algorithm failed to converge after 20 iterations, due to a signal saturation caused by strong absorption at 317 nm. In these relatively rare cases, it is beneficial to use longer wavelengths (e.g., > 320 nm) for SO 2 retrievals Theys et al., 2015), which are available from the current hyperspectral instruments such as OMI and OMPS but not TOMS.
We also calculated the residual at 312 nm (res 312 = N m − N c ), defined here as the difference between the "measured" synthetic N m and the N c at 312 nm using MS_SO2 retrieved ozone and SO 2 column amounts. Note that the 312 nm channel was not used in the MS_SO2 algorithm, and the residuals at other wavelengths are essentially zero since we are retrieving four parameters from four wavelengths. As shown in Tables S4-S7, a positive bias in the SO 2 height (when the is CMA too high when compared with the actual profile) Figure 8. Comparison of the OMPS retrieval SO 2 against the GEOS-5 model SO 2 . The TOA radiances for the OMPS retrieval were generated assuming no aerosol (a), only sulfate aerosols (b), and both ash and sulfate aerosols (c). leads to negative residuals at 312 nm, whereas a negative bias in a priori profile (CMA too low) causes positive residuals. The residuals are generally within 1-2 N value (2 %-5 % error in radiance) for SO 2 column amounts of 50-100 DU but can reach 3-7 N value (6 %-15 %) for large SO 2 amounts of 200-300 DU. While the 312 nm channel may potentially be used to retrieve SO 2 plume height for large volcanic eruptions, it is strongly affected by volcanic aerosols as demonstrated in the next section.

Ash and sulfate aerosol effects on MS_SO2 retrievals
To test the sensitivity of the MS_SO2 algorithm to ash and sulfate aerosols, an Observing System Simulation Experiment (OSSE) was conducted. The experiment used the GEOS-5 Earth system model (Molod et al., 2012;Buchard et al., 2017;Colarco et al., 2012), coupled with online Goddard Chemistry Aerosol and Radiation (GOCART) (Chin et al., 2000;Colarco et al., 2010) and Community Aerosol and Radiation Model for Atmospheres (CARMA) (Toon et al., 1988;Colarco et al., 2014). In this experiment, we considered three separate cases for a Pinatubo-like eruption scenario: (1) 12 Mt of SO 2 and no aerosols; (2) 12 Mt of SO 2 and 4 Mt of sulfate aerosols (as reported by Guo et al., 2004); and (3) 12 Mt of SO 2 , 4 Mt of sulfate aerosols, and 5 Mt of ash uniformly distributed between 18 and 22 km above the location of Pinatubo volcano, on 15 June 1991, from 06:00 to 15:00 UTC. The GEOS-5 simulated 4-D profiles of ozone, SO 2 , sulfate aerosols and volcanic ash were used as input to a VLIDORT RT model (Spurr, 2008). The model generated synthetic radiances at 317, 331, 340 and 380 nm TOMS bands, using the actual Suomi National Polar Partnership (SNPP)/OMPS-NM viewing geometry, assuming cloud-free conditions. The synthetic radiances produced by the VLIDORT were used as input to the MS_SO2 algorithm to generate "retrieved" columns of ozone and SO 2 . We note that MS_SO2 algorithm uses LUTs produced using a different TOMRAD RT model. Figure 8 compares retrieved versus true SO 2 column amounts for the three cases considered. The retrieval bias is inferred from the differences between the model SO 2 input and the SO 2 retrieved by MS_SO2, using the radiances from the model run. The no aerosol case confirms unbiased SO 2 retrievals for SO 2 column amounts less than ∼ 150 DU and small positive bias for larger SO 2 amounts. For aerosol cases where sulfates and ash were included in the simulation, we observe a negative bias for SO 2 column amounts exceeding ∼ 100 DU. These negative biases (retrieval saturation) are expected as the MS_SO2 forward model does not explicitly account for volcanic aerosols. This OSSE experiment shows the effects of heavy aerosol loading on the retrieval but also increases confidence in MS_SO2 retrievals between 15 and 100 DU, under nominal conditions, even in the presence of high aerosol concentrations.

Comparison with PCA SO 2 retrievals
We directly compared MS_SO2 retrievals with the principal component analysis (PCA) SO 2 algorithm adapted to the TOMS 6 spectral channels. In the PCA approach (Li et al., 2013, a set of principal components (PCs) is first extracted from the measured radiances using a PCA technique and ranked in descending order according to the spectral variance they each explain. If derived from SO 2 -free areas, these PCs represent geophysical processes (e.g., ozone absorption) and measurement details (e.g., wavelength shift) that are unrelated to SO 2 but may interfere with SO 2 retrievals. Next, we fit the first n ν (non-SO 2 ) PCs and the SO 2 Jacobians (∂N/∂ ) to the measured radiances (in N value) described in Eq. (10). This allows us to simultaneously estimate the coefficients of the PCs (ω) and SO 2 column amount and helps to minimize the impacts of various interfering processes: A more detailed introduction to the PCA SO 2 retrieval technique for hyperspectral instruments such as the Ozone Monitoring Instrument (OMI) and the Ozone Mapping and Profiler Suite Nadir Mapper (OMPS-NM) can be found elsewhere (e.g., Li et al., 2013Li et al., , 2017Zhang et al., 2017).
For this comparison we adapt the PCA to the discrete wavelength of N7/TOMS. The Nimbus-7 TOMS PCA SO 2 algorithm is similar to the OMI and OMPS-NM version in terms of its overall structure but differs in some implementation details. Specifically, unlike the OMI/OMPS volcanic SO 2 retrievals that use a dynamic spectral fitting window , the TOMS PCA SO 2 algorithm uses all six wavelengths available from TOMS in fitting. Also due to the small number of wavelengths, in the TOMS PCA SO 2 algorithm, we always use n ν = 5 PCs in Eq. (10), lower than the number of PCs used for OMI (n ν ≤ 20) or OMPS (n ν ≤ 15). For OMI and OMPS retrievals, SLER is derived at three wavelengths (342, 354 and 367 nm) and extrapolated to other wavelengths using a second-degree polynomial function fitted to these three wavelengths. As for TOMS, SLER is determined at 340 and 380 nm and extrapolated linearly. Additionally, while the Jacobian lookup tables are constructed using the VLIDORT radiative transfer code (Spurr, 2008) for both OMI/OMPS and Nimbus-7 TOMS, different, instrument-specific slit functions are used to band-pass the SO 2 Jacobians from the lookup tables.
We compared retrievals from the two algorithms for the first 6 d of Mount Pinatubo eruption (16-21 June 1991). The Pinatubo case provides a large sample of FoVs spanning a broad range of SO 2 amounts from 15 DU (minimum threshold) to over 400 DU. In this test of the algorithm, MS_SO2 and PCA retrievals were generated assuming a CMA of 18 km.

June 1991 eruption of Mount Pinatubo
Mount Pinatubo is a large stratovolcano located at 15 • 08 N, 120 • 21 E in western Luzon, Philippines, that erupted explosively on 15 June 1991, following weeks of precursory activity. TOMS SO 2 imagery on 15 June shows a narrow, elongated SO 2 ash plume extending to the west from the location of the volcano. On the following day TOMS measured a massive SO 2 plume to the west of the volcano (Bluth et al., 1992). TOMS continued tracking the daily evolution of the Pinatubo volcanic cloud as it encircled the Earth over a period of about 22 d. Previous estimates of the Pinatubo SO 2 height (CMA) range between 18 and 25 km (Self et al., 1996;Guo et al., 2004). Figure 9 shows TOMS daily SO 2 maps produced with the MS_SO2 and the PCA algorithms for the 6 d period from 16 to 21 June. Corresponding ash index (AI) imagery from MS_SO2 is shown in Fig. 10. SO 2 and AI imagery for 16 June show a large SO 2 ash cloud propagating to the west. AI values range from 4 to above 12 across the plume. The AI values decreased over the following days due to wind advection and wet deposition (Guo et al., 2004). As the SO 2 cloud area continues to expand, total SO 2 mass remains high, while SO 2 peak values decrease, which is expected from cloud dispersion. The MS_SO2 and the PCA imagery show excellent qualitative agreement in resolving the plume area and internal SO 2 plume structure, as inferred from the SO 2 gradients across the peak regions of the cloud. Note that for 16-19 June, part of the observed cloud is missing due to a known mechanical problem with the TOMS instrument. These missing regions can be clearly identified in the imagery. Figure 11 shows a scatterplot comparing the MS_SO2 and PCA retrievals for the 6 d time series, which included over 7000 matching FoVs. These results show the retrievals are in close quantitative agreement, with a correlation of 0.993 and a slope of 1.00. Since the two algorithms apply fundamentally different approaches to retrieving SO 2 , this level of agreement is impressive considered over such a broad range of values.
We further compared quantitative estimates of SO 2 cloud mass, peak SO 2 and plume area. For this comparison, we also considered results from the Krueger-Kerr algorithm (KK), based on the published results of Guo et al. (2004). Table 3 displays daily estimates of the SO 2 cloud mass and peak SO 2 amounts for the MS_SO2, PCA and KK algorithms for the 6 d period. Guo et al. (2004) applied a modified version of the KK algorithm that assumes a radiative transfer air mass factor (AMF), which accounts for the a priori ozone and SO 2 absorption profiles (Krotkov et al., 1997). The early SO 2 mass estimates by Bluth et al. (1992) derived from Pinatubo eruption assumed a geometrical AMF. Also note that Guo et al. (2004) interpolated across the missing data regions of the plume on 16, 18 and 19 June using a punctual kriging statistical analysis. Here, we did not correct for the missing data. The three algorithms are in good overall agreement for the period from 17 to 21 June, with the differences within 10 % compared to MS_SO2. The most significant differences between the three algorithms are observed on 16 June under conditions of heavy ash loading. KK mass tonnage estimates exceeded MS_SO2 by over 24 % even though MS_SO2 and the PCA differ by just 2 %. Some of the difference between KK and the other two algorithms can be attributed to the fact that the Guo et al. (2004) estimates include contribution from the missing data region at the northern boundary of the plume (compare SO 2 and aerosol imagery), but this contribution does not nearly account for the total difference in Table 3.
The differences can be explained by considering how each algorithm is affected by aerosols. MS_SO2 accounts for ash by retrieving the spectral dependence at 340 nm, which is then adjusted iteratively to correct the reflectivity at the two absorbing channels. As explained in Sect. 3.2, absorbing aerosols in the column can cause possible ozone anomalies, which decrease . The KK algorithm (Krueger et al., 1995) accounts for ash implicitly by retrieving two linear spectral parameters that adjust calculated N c to match measured N m . Like MS_SO2, the KK radiative path LUTs are based on TOMRAD calculations that do not explicitly account for ash (Krotkov et al., 1997). Krueger et al. (1995) estimated that ash aerosols can cause errors in the retrieval up to +30%, depending on the ash size distribution. The PCA algorithm,    in contrast, accounts for ash in the separation and ordering of the principal components. The differences between MS_SO2 and KK on 16 and 17 June can be partly ascribed to the effects of aerosols on the retrievals. By 18 June, the ash and SO 2 clouds have mostly separated, though, aerosol indices over 4 are still observed in some regions of the plume. Pinatubo did not erupt again after the major eruption on 15 June, yet the three algorithms show retrieved SO 2 mass increases on 17 and 20 June (the PCA and KK retrievals also indicate a small increase on 18 June). Guo et al. (2004) attribute these increases to the sequestering of volcanic SO 2 by ice-ash mixtures in the plume. They propose the sequestered SO 2 was released at a later time through sublimation of ice in the lower stratosphere. The oxidation of hydrogen sulfide offers another mechanism to account for the observed mass increases in the days following the eruption. The combined results of the three algorithms support the conclusions of Guo et al. (2004) that the observed mass increases in the temporal evolution of the plume are real.
Overall, the PCA retrieved 3 % more total mass tonnage than MS_SO2. These differences are attributed to differences in how the MS_SO2 algorithm handles aerosols and differences in the area of the plume due to differences in the retrieval near the sensitivity threshold (∼ 15 DU). Ash, sulfates and high SO 2 amounts impact the ozone retrieval, for as was seen in Sect. 3.2, systematic errors in SO 2 are anticorrelated with errors in O 3 (see Fig. S1). For the case of the KK algorithm, the total ozone retrieved inside the SO 2 plume can be unrealistically low and even negative in an extreme event like Mount Pinatubo shown in Figs. S4 and S5. Figure S4 compares the KK ozone retrieval with MS_SO2 step 2 ozone retrieval and Fig. S5 compares scatterplots of SO 2 and total ozone for 17 and 18 June. Table 4 provides estimates of the plume area for the MS_SO2 and PCA. The area of the plume is most sensitive to the minimum detection threshold around the edges of the SO 2 cloud. MS_SO2 and the PCA algorithms were directly compared by computing the areal sum of all the pixels where > 15 DU (Fig. 9). For the 6 d study period, the plume increased in size from about a little over 2 × 10 6 km 2 to ∼ 9 × 10 6 km 2 The PCA trends observed a larger cloud area for five of the 6 d, with most of the observed differences within 7 %. On 16 June, shortly after the major eruption of 15 June, the estimated area for the PCA is about 15 % greater than for MS_SO2. The fresh plumes are opaque, which result in underestimating of SO 2 mass by all BUV algorithms due to the mixing of aerosols (Krotkov et al., 1997). The PCA appears slightly more sensitive to SO 2 near the edges of the cloud, where aerosol loading is high (AI > 1.5). It should be noted that the soft calibration applied to the 340 nm channel, described in Sects. 3.3 and S3.3, may also contribute to the lowering the sensitivity around the edges of the plume. This correction effectively lowered the background by about 3 DU.

Conclusions
This paper describes a discrete multi-satellite UV wavelength algorithm (MS_SO2) for retrieving volcanic SO 2 that was used operationally to process measurements from the heritage Nimbus-7 TOMS and the Deep Space Climate Observatory Earth Polychromatic Imaging Camera Marshak et al., 2018). The MS_SO2 algorithm retrieves four parameters (SO 2 , O 3 , dR/dλ and R s ) and can be used to process data from current hyperspectral UV spectrometers, such as SNPP/OMPS and Aura/OMI, using a convolved, discrete set of wavelengths, offering a viable means for intercomparing volcanic SO 2 retrievals from different missions. We estimated random (noise) and systematic errors, related to the effects of volcanic aerosols and uncertainties in SO 2 height and partly corrected for absorbing ash, using positive aerosol index (AI) as a proxy for applying a Step 2 correction to the SO 2 retrievals. The correction could still underestimate SO 2 mass during the first days after extremely large eruptions (VEI > 3) due to BUV saturation. In such cases we recommend estimating e-folding time of the SO 2 decay, using later measurements and extrapolating SO 2 mass exponentially back in time to the eruption day .
The TOMS Observing System Simulation Experiment simulation, using synthetic radiances, shows unbiased MS_SO2 retrievals of for SO 2 < 100-150 DU but low biases for larger SO 2 amounts due to the presence of ash and sulfate aerosols. Therefore, operational MS_SO2 retrievals should provide a low boundary constraint on the SO 2 mass injected into the atmosphere from large eruptions during first days after an eruption. The algorithm can be further improved by explicitly accounting for volcanic ash and sulfate aerosols, which was not feasible in the operational processing.
The MS_SO2 retrieval is also sensitive to differences between the a priori and actual SO 2 center of mass altitude. Since this key parameter is not retrieved, the TOMS SO 2 product provides separate SO 2 column amounts assuming three different SO 2 altitudes (8, 13 and 18 km). Users should base their analysis on the altitude that is most appropriate for a particular eruption.
To assess the overall accuracy of the TOMS SO 2 retrievals, we compared MS_SO2 and independent PCA algorithms for the first 6 d following the 1991 Pinatubo eruption. The daily time series of SO 2 retrievals showed high correlation (R 2 = 0.986) and excellent agreement between the two retrievals over a broad SO 2 range between 15 and 400 DU. We also compared the SO 2 mass, peak SO 2 amounts and plume area with the heritage Krueger-Kerr algorithm. This threeway comparison showed the SO 2 mass within 10 % for all days, except on 16 June, when the Krueger-Kerr algorithm retrieved 24 % higher SO 2 mass. This could be explained by interpolation over a region of missing TOMS measurements on 16 June (Guo et al., 2004). The remaining differences between current MS_SO2 and the PCA algorithms (3 %-7 %) are attributed to the differences in handling of aerosols and different sensitivity thresholds of the algorithms.
The reprocessed Nimbus-7 TOMS volcanic SO 2 data set (TOMSN7SO2) is now publicly available through the Goddard Earth Sciences Data and Information Services Center (GES DISC) as part of the NASA's Making Earth System Data Records for Use in Research Environments (MEa-SUREs) program (Krotkov et al., 2019). We plan to reprocess all follow-up multispectral UV (TOMS) and hyperspectral UV (OMI, OMPS) missions ( Fig. 1) with MS_SO2 and PCA algorithms to keep updating our multi-satellite volcanic SO 2 mass database archived at GES DISC (Carn, 2019). It is important to continue quantifying SO 2 emissions from small explosive eruptions, as they may, collectively, play an important role in sustaining the persistent, background stratospheric aerosol layer, which is an important factor in global climate forcing. Data availability. Our data can be publicly accessed at https://disc.gsfc.nasa.gov/datasets/TOMSN7SO2_3/summary? keywords=TOMSSO2 (Krotkov et al., 2019).
Author contributions. BF implemented the MS_SO2 algorithm, processed and archived the N7TOMS SO 2 data record, wrote most of the text and supplement, produced the figures, and coordinated with the other co-authors. NK advised on paper organization and figures and wrote Sect. 2.2 and the Conclusion. PB developed the theoretical basis for the MS_SO2 algorithm and wrote Sect. 2.1 on the heritage ozone algorithm. CL provided PCA SO 2 results and contributed to Sects. 4.2 and 5 and the Supplement. SC wrote the Introduction and important scientific analysis of the data. EH contributed to Sect. 4.2.2, quantifying the effects of ash and sulfate aerosols on the MS_SO2 retrieval. PL proposed and implemented the data format used in the MEaSUREs archived SO 2 products.
Competing interests. The authors declare that they have no conflict of interest.
Financial support. This work was supported by NASA's Making Earth System data records for Use in Research Environments (NNH12ZDA001N-MEASURES) program. Can Li acknowledges support from NASA Earth Science Division for development and analysis of hyperspectral PCA SO 2 products for OMPS (grant no. 80NSSC18K0688). Eric Hughes was supported by NASA grant to the University of Maryland no. NNX13AG51G. Simon A. Carn acknowledges support from NASA grant NNX13AF50G.
Review statement. This paper was edited by Michel Van Roozendael and reviewed by two anonymous referees.