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

This paper describes a new discrete wavelength algorithm developed for retrieving volcanic sulfur dioxide (SO₂) vertical column density (VCD) from UV observing satellites. The Multi-Satellite SO₂ 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 SO₂ 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 SO₂ 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 SO₂ VCDs from contemporary hyperspectral UV instruments, such as OMI and OMPS, to create consistent, multi-satellite, long-term volcanic SO₂ climate data records.

1 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₂) 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₂ 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₂ emissions are essential.

The satellite record of volcanic SO₂ 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., 2003, 2016; Carn, 2019; Prata et al., 2003). Measurements in the UV band have a longer heritage, since the first satellite detection of volcanic SO₂ 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₂ 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₂ column than IR retrievals, and hence the former have been the mainstay of volcanic SO₂ monitoring during the satellite era to date. The volcanic SO₂ climatology from 1978 to the present (Fig. 1, Carn, 2019) reveals highly variable inter-annual volcanic SO₂ forcing dominated by two major eruptions (El Chichón in 1982 and Pinatubo 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).

Figure 1 Multi-decadal record of SO₂ 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₂ production (omitting SO₂ 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).

One of the key challenges in assembling a long-term, consistent, satellite-based volcanic SO₂ 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₂ 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₂ sensitivity occurred when the multispectral, six-channel TOMS instruments were superseded by hyperspectral UV sensors, such as the Global Ozone Monitoring Experiment (GOME, 1995–2003; Khokhar et al., 2005), the Scanning Imaging Absorption Spectrometer for Atmospheric Chartography (SCIAMACHY, 2002–2012; Lee et al., 2008), the Ozone Monitoring Instrument (OMI, 2004–ongoing; Krotkov et al., 2006), the Ozone Mapping and Profiler Suite (OMPS, 2012–ongoing; 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₂ 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₂ 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., 1995, 2000; Krotkov et al., 1997, 2006; Yang et al., 2007, 2010; Li et al., 2013, 2017; Theys 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₂ climatology across multiple UV satellite missions. In this paper we describe a new multi-satellite SO₂ 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₂ record, we apply the MS_SO2 algorithm to the Nimbus-7 TOMS (N7/TOMS) measurements (1978–1993) and present a reanalysis of some of the most significant eruptions of the N7/TOMS mission.

2 Heritage satellite ozone and SO₂ algorithms

Ozone and SO₂ 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₃ and SO₂ cross sections and the SO₂∕O₃ 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₂ molecule is 2.5 times more absorbing than O₃ at 317 nm, while O₃ is 6 times more absorbing at 331 nm. These differences allow for simultaneous multispectral retrievals of O₃ and SO₂.

Figure 2 Spectral dependence of laboratory-measured SO₂ (black) and O₃ (red) cross sections between 310 and 340 nm at TOMS FWHM ∼1 nm. The SO₂∕O₃ ratio (green) is shown with the scale on the right axis. The nominal locations of the N7/TOMS absorbing bands (317, 331, 340 nm) are shown by vertical blue lines (blue).

2.1 Heritage BUV ozone algorithms

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₂ 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

$$\begin{array}{}\text{(1)}& N=-\mathrm{100}{\log }_{\mathrm{10}}\left({\displaystyle \frac{I}{F}}\right).\end{array}$$

F is the extraterrestrial solar irradiance. By linearly interpolating N between total ozone nodes, one forms the N-Ω curves 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₂ 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 increased 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₂ 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).

2.2 Heritage TOMS SO₂ algorithms

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₂ released into the atmosphere during the event. To estimate the SO₂ inside the plume region, he separated the SO₂ and O₃ 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₂ 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₃ and SO₂ radiance contributions, based on an earlier methodology developed by Kerr et al. (1980) to retrieve the SO₂ column from the ground with a Brewer spectrophotometer. This method assumed that the BUV radiation was attenuated by the two absorbing species (O₃, SO₂), leading to an equation describing BUV radiance, I, for a given wavelength, λ, corresponding to the TOMS field of view (FoV):

$$\begin{array}{}\text{(2)}& I\left(\mathit{\lambda }\right)=aF\left(\mathit{\lambda }\right)\exp \left[-b\mathit{\lambda }+S_{\mathrm{g}}({\mathit{\tau }}_{{\mathrm{O}}_{\mathrm{3}}}+{\mathit{\tau }}_{{\mathrm{SO}}_{\mathrm{2}}})\right].\end{array}$$

In Eq. (2), F is the incoming solar flux, S_g is the geometrical optical path (air mass factor, AMF), and ${\mathit{\tau }}_{{\mathrm{O}}_{\mathrm{3}}}$ and ${\mathit{\tau }}_{{\mathrm{SO}}_{\mathrm{2}}}$ are the vertical optical thicknesses for O₃ and SO₂, 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₂ and O₃ vertical column 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₂ 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₂ mass tonnage from medium to large explosive eruptions using TOMS BUV measurements (Krueger et al., 2000; Carn et al., 2003).

3 New MS_SO2 algorithm

The new discrete wavelength SO₂ algorithm (MS_SO2) builds on the heritage of the TOMS total ozone algorithm (Sect. 2.1) but adds sulfur dioxide (SO₂) 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₂ 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). TOMRAD accounts for all orders of polarized Rayleigh scattering and for the gaseous absorption (e.g., O₃ and SO₂), 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.

Table 1 Input parameters used in construction of the Nimbus-7 TOMS LUTs.

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, N7TOMS-specific 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₃ and SO₂. The LUTs include 21 total ozone nodes and 12 total SO₂ nodes for each of the three assumed SO₂ 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₂ concentration. LUTs for SO₂ are generated for three different CMAs: 8 km (middle troposphere, TRM), 13 km (upper tropospheric, TRU) and 18 km (lower stratospheric, STL). A constant standard deviation of σ=2 km is assumed for each SO₂ profile.

The MS_SO2 algorithm retrieves a four-parameter state vector, x,defined below as

$$\begin{array}{}\text{(3)}& \mathit{x}=\left(\begin{array}{c}\mathrm{\Sigma }\\ \begin{array}{c}\mathrm{\Omega }\\ \mathrm{d}Rs/\mathrm{d}\mathit{\lambda }\\ R_{\mathrm{s}}\end{array}\end{array}\right),\end{array}$$

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.

3.1 Step 1 retrieval

Our step 1 inversion starts with an initial state vector x₀, consisting of first guesses for Σ₀, Ω₀ 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:

$$\begin{array}{}\text{(4)}& \mathrm{d}\mathbf{N}=\mathbf{K}\mathrm{d}\mathit{x},\end{array}$$

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:

$$\begin{array}{}\text{(5)}& K_{i,j}={\displaystyle \frac{\partial N_{\mathrm{c},i}}{\partial x_j}},i,j=\mathrm{1},\mathrm{4},\end{array}$$

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₃ and SO₂ cross sections are negligible at 380 nm, the R_s and $\partial N_{\mathrm{380}}/\partial R_{\mathrm{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, ${\mathit{x}}_k={\mathit{x}}_{k-\mathrm{1}}+\text{d}{\mathit{x}}_k$, $k=\mathrm{1},\mathrm{2}$,… until it converges on a solution as described below:

$$\begin{array}{}\text{(6a)}& {\displaystyle }\mathrm{d}{\mathbf{N}}_k={\mathbf{N}}_m-{\mathbf{N}}_{c,k-\mathrm{1}}={\mathbf{K}}_{k-\mathrm{1}}\mathrm{d}{\mathit{x}}_k,\text{(6b)}& {\displaystyle }\mathrm{d}{\mathit{x}}_k={\mathit{x}}_k-{\mathit{x}}_{k-\mathrm{1}}={\mathbf{K}}_{k-\mathrm{1}}^{-\mathrm{1}}\mathrm{d}{\mathbf{N}}_k.\end{array}$$

Since O₃ and SO₂ exhibit small absorption at 340 nm, a nonzero R spectral slope (i.e., $\mathrm{d}R/\mathrm{d}\mathit{\lambda }\ne \mathrm{0}$) accounts for the radiative effects of aerosols and surface reflectance (e.g., sun glint).

Table 2 Retrieved state vector. n/a – not applicable.

^*Ω₀ is a climatological value for each of the three latitude bands.

As indicated in Table 2, the algorithm initially assumes zero R-λ dependence (i.e., $\mathrm{d}R/\mathrm{d}\mathit{\lambda }=\mathrm{0}$); however, absorbing aerosols (smoke, dust and volcanic ash) cause $\mathrm{d}R/\mathrm{d}\mathit{\lambda }\ne \mathrm{0}$.

The algorithm uses retrieved spectral slope dR∕dλ in Eq. (7) below to update the calculated LERs after each iteration:

$$\begin{array}{}\text{(7)}& R_j=R_{\mathrm{s}}+{\displaystyle \frac{\partial R}{\partial \mathit{\lambda }}}\times \left({\mathit{\lambda }}_j-{\mathit{\lambda }}_R\right)j=\mathrm{1},\mathrm{4},\end{array}$$

where λ_j=312, 317, 331, 340 nm and λ_R=380 nm. When SO₂ 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₃₁₂ 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.

3.2 Step 2 retrieval

The MS_SO2 forward model accounts for O₃ and SO₂ 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₂ 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₂ 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:

$$\begin{array}{}\text{(8)}& \mathrm{AI}={\displaystyle \frac{\partial N_{\mathrm{340}}}{\partial R}}{\displaystyle \frac{\partial R}{\partial \mathit{\lambda }}}\left({\mathit{\lambda }}_{\mathrm{340}}-{\mathit{\lambda }}_{\mathrm{380}}\right)=-\mathrm{40}\cdot {\displaystyle \frac{\partial N_{\mathrm{340}}}{\partial R}}{\displaystyle \frac{\partial R}{\partial \mathit{\lambda }}}.\end{array}$$

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₂ > 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₂ can be less than 15 DU due to high aerosol concentrations. In this case, it is assumed that the step 1 SO₂ may have been underestimated due to the ozone error caused by high aerosol concentrations (in these cases, the SO₂ 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₂ 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₂ 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₂ (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₂ mass. In this case, the SO₂ 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.

Figure 3 Aerosol index for the El Chichón eruption on 4 April 1982, computed from retrieved dR∕dλ.

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₂ cloud mass is insignificant.

Figure 4 MS_SO2 maps showing (a) step 1 total column O₃, (b) step 2 total column O₃, (c) step 1 total column SO₂ and (d) step 2 total column SO₂ from the El Chichón eruption on 4 April 1982.

Step 2 follows a methodology similar to the original residual method developed by Krueger (1983), which separated the O₃ and SO₂ contributions by subtracting the measured BUV reflectance in the unperturbed region from the BUV radiance anomaly associated with the SO₂ 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₂ 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₂ 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₂ to sulfate conversion e-folding lifetime using weeks of measurements of the total SO₂ cloud daily mass and extrapolating it back in time to estimate total SO₂ mass emitted on eruption day. This “day one” time extrapolated SO₂ mass is typically larger than retrieved on days immediately following the eruption (Krotkov et al., 2010).

3.3 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₂ 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₂, Σ_true and the mean of the distribution, 〈Σ〉_clean, to equal zero such that

$$\begin{array}{}\text{(9)}& {\mathrm{\Sigma }}_{\mathrm{true}}=\langle \mathrm{\Sigma }{\rangle }_{\mathrm{clean}}=\mathrm{0}.\end{array}$$

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₃₄₀ 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₂ before (dashed) and after (solid) applying the correction for 11 November 1981. The mean bias is reduced to < 1 DU after applying the correction.

Figure 5 Probability density function of SO₂ background before (dashed) and after (solid) applying N₃₄₀ value correction.

4 Error analysis

4.1 Random errors and SO₂ detection limit

The random errors in the MS_SO2 retrieval were estimated from the standard deviation in the SO₂ 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 deviations 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₂ 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 Polar-orbiting 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 6 (a) PDF of SO₂ background (noise distribution) for TOMS and OMPS based on orbits from clean regions of the central Pacific, and (b) standard deviations of background SO₂ for TOMS and the OMPS nadir mapper as a function of the swath position. OMPS noise is more than a factor of 2–3 lower than TOMS and less dependent on cross-track position.

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₂ 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⁻¹) (Zhang et al., 2017).

4.2 Systematic errors in volcanic SO₂ plumes

In this section, we evaluate systematic errors of the MS_SO2 retrievals of volcanic SO₂. The two most significant errors are caused by volcanic aerosols (ash and sulfate) and incorrect assumptions regarding the SO₂ profile, namely the plume height. The radiance tables used by the algorithm account for ozone and SO₂ 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 (Rose et al., 2003; Guo et al., 2004). The precomputed radiance tables used by MS_SO2 assume an SO₂ column amount and an a priori CMA and standard deviation (Sect. 3). An incorrect CMA assumption can cause significant SO₂ errors that vary with viewing geometry, ozone and SO₂ column amounts. We characterize these error sources by applying the MS_SO2 algorithm to synthetic radiances.

Figure 7 VLIDORT calculated SO₂ column Jacobians ($\partial N/\partial \mathrm{\Sigma }$) at 317 nm for typical conditions in the tropics (SZA $=\mathrm{10}{}^{\circ }$, relative azimuthal angle, RAA $=\mathrm{90}{}^{\circ }$, O₃=275 DU) but different SO₂ column amounts (50, 100, 200 and 300 DU) and center mass altitudes (11–20 km). For these calculations, Gaussian SO₂ profiles with the same standard deviation (2 km) were assumed: (a) VZA =0 and R=0.05, (b) VZA =0 and R=0.50, (c) VZA =60 and R=0.05, and (d) VZA =60 and R=0.50.

4.2.1 Uncertainties due to SO₂ plume height

To understand retrieval errors in MS_SO2 algorithm due to assumed a priori SO₂ 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₂ Jacobians $\partial N/\partial \mathrm{\Sigma }$ at 317 nm for different SO₂ amounts, Σ, nadir angles and scene reflectance as function of the assumed SO₂ height (center of mass altitude, CMA). The Jacobians generally increase with the CMA, meaning that satellite BUV measurements are more sensitive to SO₂ at higher altitudes. This means that the MS_SO2 algorithm will overestimate (underestimate) the SO₂ column amount if the CMA of the a priori profile is lower (higher) than that of the actual SO₂ profile. On the other hand, the sensitivity of SO₂ Jacobians with respect to CMA is affected by several other factors, particularly SO₂ 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₂ Jacobians to CMA is greater for SO₂ plumes with large SO₂ 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₂ Jacobians to CMA. From these calculations, we can estimate the errors in the SO₂ 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₂ plumes from relatively moderate eruptions (∼50 DU), the relative errors in SO₂ Jacobians due to the error in the CMA are mostly within ±10 %. But for plumes with large SO₂ loading (∼200–300 DU) from explosive eruptions such as Pinatubo, the relative error in SO₂ 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₂ 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₂ and O₃ column amounts were compared with assumed in VLIDORT calculations (Tables S4–S7). As shown in the tables, for SO₂ plumes with a modest loading (∼50 DU), the relative errors in SO₂ column amounts, due to a 2 km error in the a priori profile are typically 10 % or less, whereas the relative errors in O₃ are within 1 %. For plumes with large SO₂ loadings (200–300 DU), the errors in SO₂ 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₂=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₂ retrievals (Li et al., 2017; Theys et al., 2015), which are available from the current hyperspectral instruments such as OMI and OMPS but not TOMS.

Figure 8 Comparison of the OMPS retrieval SO₂ against the GEOS-5 model SO₂. The TOA radiances for the OMPS retrieval were generated assuming no aerosol (a), only sulfate aerosols (b), and both ash and sulfate aerosols (c).

We also calculated the residual at 312 nm (res${}_{\mathrm{312}}=N_{\mathrm{m}}-N_{\mathrm{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₂ 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₂ height (when the is CMA too high when compared with the actual profile) 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₂ column amounts of 50–100 DU but can reach 3–7 N value (6 %–15 %) for large SO₂ amounts of 200–300 DU. While the 312 nm channel may potentially be used to retrieve SO₂ plume height for large volcanic eruptions, it is strongly affected by volcanic aerosols as demonstrated in the next section.

4.2.2 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₂ and no aerosols; (2) 12 Mt of SO₂ and 4 Mt of sulfate aerosols (as reported by Guo et al., 2004); and (3) 12 Mt of SO₂, 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₂, 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₂. We note that MS_SO2 algorithm uses LUTs produced using a different TOMRAD RT model.

Figure 8 compares retrieved versus true SO₂ column amounts for the three cases considered. The retrieval bias is inferred from the differences between the model SO₂ input and the SO₂ retrieved by MS_SO2, using the radiances from the model run. The no aerosol case confirms unbiased SO₂ retrievals for SO₂ column amounts less than ∼150 DU and small positive bias for larger SO₂ amounts. For aerosol cases where sulfates and ash were included in the simulation, we observe a negative bias for SO₂ 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.

5 Comparison with PCA SO₂ retrievals

We directly compared MS_SO2 retrievals with the principal component analysis (PCA) SO₂ algorithm adapted to the TOMS 6 spectral channels. In the PCA approach (Li et al., 2013, 2017), 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₂-free areas, these PCs represent geophysical processes (e.g., ozone absorption) and measurement details (e.g., wavelength shift) that are unrelated to SO₂ but may interfere with SO₂ retrievals. Next, we fit the first n_ν (non-SO₂) PCs and the SO₂ Jacobians ($\partial N/\partial \mathrm{\Sigma }$) to the measured radiances (in N value) described in Eq. (10). This allows us to simultaneously estimate the coefficients of the PCs (ω) and SO₂ column amount and helps to minimize the impacts of various interfering processes:

$$\begin{array}{}\text{(10)}& N\left(\mathit{\omega },\mathrm{\Sigma }\right)={\sum }_{i=\mathrm{1}}^{n_v}{\mathit{\omega }}_i{v}_i+{\displaystyle \frac{\partial N}{\partial \mathrm{\Sigma }}}\mathrm{\Sigma }.\end{array}$$

A more detailed introduction to the PCA SO₂ 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., 2013, 2017; Zhang et al., 2017).

For this comparison we adapt the PCA to the discrete wavelength of N7/TOMS. The Nimbus-7 TOMS PCA SO₂ 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₂ retrievals that use a dynamic spectral fitting window (Li et al., 2017), the TOMS PCA SO₂ algorithm uses all six wavelengths available from TOMS in fitting. Also due to the small number of wavelengths, in the TOMS PCA SO₂ 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₂ 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₂ 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₂ imagery on 15 June shows a narrow, elongated SO₂ ash plume extending to the west from the location of the volcano. On the following day TOMS measured a massive SO₂ 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₂ height (CMA) range between 18 and 25 km (Self et al., 1996; Guo et al., 2004).

Figure 9 Daily SO₂ imagery for MS_SO2 and the PCA using data from TOMS overpasses of the Pinatubo eruption cloud between 15 and 21 June 1991: (a) MS_SO2 for 16 June and (b) MS_SO2 for 17 June, (c) MS_SO2 for 18 June and (d) PCA for 16 June, (e) PCA for 17 June and (f) PCA for 18 June, (g) MS_SO2 for 19 June and (h) MS_SO2 for 20 June, (i) MS_SO2 for 21 June and (j) PCA for 19 June, and (k) PCA for 20 June and (l) PCA for 21 June.

Figure 9 shows TOMS daily SO₂ 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₂ and AI imagery for 16 June show a large SO₂ 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₂ cloud area continues to expand, total SO₂ mass remains high, while SO₂ 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₂ plume structure, as inferred from the SO₂ 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 10 Daily AI imagery retrieved using MS_SO2 between 16 and 21 June 1991. Contours show SO₂ levels from Fig. 9. Positive AI values over India and the Arabian peninsula are due to dust aerosols, not related to the Pinatubo ash cloud: (a) 16 June, (b) 17 June, (c) 18 June, (d) 19 June, (e) 20 June and (f) 21 June.

Figure 11 Scatterplot of retrieved SO₂ using PCA and MS_SO2 algorithms for the period 16–21 June 1991.

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₂, this level of agreement is impressive considered over such a broad range of values.

Table 3 Daily SO₂ mass and maximal SO₂ values for MS_SO2, PCA and KK algorithms for the 6 d following the Pinatubo eruption on 15 June 1991.

^* Guo et al. (2004) interpolated values in the missing data region seen in maps for 16, 18 and 19 June.

We further compared quantitative estimates of SO₂ cloud mass, peak SO₂ 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₂ cloud mass and peak SO₂ 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₂ absorption profiles (Krotkov et al., 1997). The early SO₂ 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₂ 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₂ 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₂ 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₂ by ice–ash mixtures in the plume. They propose the sequestered SO₂ 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₂ amounts impact the ozone retrieval, for as was seen in Sect. 3.2, systematic errors in SO₂ are anticorrelated with errors in O₃ (see Fig. S1). For the case of the KK algorithm, the total ozone retrieved inside the SO₂ 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₂ and total ozone for 17 and 18 June.

Table 4 SO₂ plume area and number of fields of view (NFoVs) where the retrieved SO₂ exceeded 15 DU using the MS_SO2 and PCA algorithms for the 6 d following for the Pinatubo eruption on 15 June 1991.

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₂ 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⁶ km² to $\sim \mathrm{9}\times {\mathrm{10}}^{\mathrm{6}}$ km² 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₂ mass by all BUV algorithms due to the mixing of aerosols (Krotkov et al., 1997). The PCA appears slightly more sensitive to SO₂ 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.

6 Conclusions

This paper describes a discrete multi-satellite UV wavelength algorithm (MS_SO2) for retrieving volcanic SO₂ that was used operationally to process measurements from the heritage Nimbus-7 TOMS and the Deep Space Climate Observatory Earth Polychromatic Imaging Camera (Carn et al., 2018; Marshak et al., 2018). The MS_SO2 algorithm retrieves four parameters (SO₂, O₃, 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₂ retrievals from different missions.

We estimated random (noise) and systematic errors, related to the effects of volcanic aerosols and uncertainties in SO₂ height and partly corrected for absorbing ash, using positive aerosol index (AI) as a proxy for applying a Step 2 correction to the SO₂ retrievals. The correction could still underestimate SO₂ 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₂ decay, using later measurements and extrapolating SO₂ mass exponentially back in time to the eruption day (Krotkov et al., 2010).

The TOMS Observing System Simulation Experiment simulation, using synthetic radiances, shows unbiased MS_SO2 retrievals of for SO₂ < 100–150 DU but low biases for larger SO₂ amounts due to the presence of ash and sulfate aerosols. Therefore, operational MS_SO2 retrievals should provide a low boundary constraint on the SO₂ 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₂ center of mass altitude. Since this key parameter is not retrieved, the TOMS SO₂ product provides separate SO₂ column amounts assuming three different SO₂ 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₂ 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₂ retrievals showed high correlation (R²=0.986) and excellent agreement between the two retrievals over a broad SO₂ range between 15 and 400 DU. We also compared the SO₂ mass, peak SO₂ amounts and plume area with the heritage Krueger–Kerr algorithm. This three-way comparison showed the SO₂ mass within 10 % for all days, except on 16 June, when the Krueger–Kerr algorithm retrieved 24 % higher SO₂ 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₂ 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 (MEaSUREs) 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₂ mass database archived at GES DISC (Carn, 2019). It is important to continue quantifying SO₂ 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.