A Novel Approach to Estimating Time-Averaged Volcanic SO₂ Fluxes from Infrared Satellite Measurements

Abstract

Long-term continuous time series of ${\mathrm{SO}}_2$ emissions are considered critical elements of both volcano monitoring and basic research into processes within magmatic systems. One highly successful framework for computing these fluxes involves reconstructing a representative time-averaged ${\mathrm{SO}}_2$ plume from which to estimate the ${\mathrm{SO}}_2$ source flux. Previous methods within this framework have used ancillary wind datasets from reanalysis or numerical weather prediction (NWP) to construct the mean plume and then again as a constrained parameter in the fitting. Additionally, traditional ${\mathrm{SO}}_2$ datasets from ultraviolet (UV) sensors lack altitude information, which must be assumed, to correctly calibrate the ${\mathrm{SO}}_2$ data and to capture the appropriate NWP wind level which can be a significant source of error. We have made novel modifications to this framework which do not rely on prior knowledge of the winds and therefore do not inherit errors associated with NWP winds. To perform the plume rotation, we modify a rudimentary computer vision algorithm designed for object detection in medical imaging to detect plume-like objects in gridded ${\mathrm{SO}}_2$ data. We then fit a solution to the general time-averaged dispersion of ${\mathrm{SO}}_2$ from a point source. We demonstrate these techniques using ${\mathrm{SO}}_2$ data generated by a newly developed probabilistic layer height and column loading algorithm designed for the Cross-track Infrared Sounder (CrIS), a hyperspectral infrared sensor aboard the Joint Polar Satellite System’s Suomi-NPP and NOAA-20 satellites. This ${\mathrm{SO}}_2$ data source is best suited to flux estimates at high-latitude volcanoes and at low-latitude, but high-altitude volcanoes. Of particular importance, IR ${\mathrm{SO}}_2$ data can fill an important data gap in the UV-based record: estimating ${\mathrm{SO}}_2$ emissions from high-latitude volcanoes through the polar winters when there is insufficient solar backscatter for UV sensors to be used.

1. Introduction

Many methods now exist to estimate volcanic sulfur dioxide (${\mathrm{SO}}_2$) fluxes from satellite data. Over the last decade, much work has focused on long term, global monitoring of continuous low–moderate level volcanic ${\mathrm{SO}}_2$ emissions. These efforts are particularly useful for monitoring trends in volcanic activity as well as constraining the global volcanic ${\mathrm{SO}}_2$ flux for weather and climate studies. These methods for estimating long-term emission rates have relied on constructing and analyzing a time-averaged volcanic plume from satellite ${\mathrm{SO}}_2$ data generated by ultraviolet (UV) sensors. To do this, either the wind field advecting ${\mathrm{SO}}_2$ emissions must be remarkably stable e.g., [1] or the many plume observations in the averaging period must be rotated into a common wind field and then averaged together e.g., [2,3,4,5]. Once the time-averaged plume is constructed, a simplified plume dispersion model is fit to the time-averaged data. As described by [2], the measurement of ${\mathrm{SO}}_2$ flux is controlled by an estimate of the cloud mass and of the decay time which are estimated from the fitting results. Since the mean cloud mass is directly accessible from satellite measurements, the decay time (or lifetime) is the primary uncertain parameter.

Typically, numerical weather prediction (NWP) or re-analysis winds are used to analyze the wind direction to determine the rotation and then the wind speeds are used in the plume model fitting process. Errors in NWP, particularly in wind direction can have a significant effect on the reconstructed plume. In this process, a plume dispersal height must be assumed at which to capture the wind direction and speed and to estimate the ${\mathrm{SO}}_2$ vertical column density (VCD) itself. Typically, individual ${\mathrm{SO}}_2$ plume observations are assumed to have been advected by straight-line winds at one altitude; however, volcanic plumes are hot-sourced and can rise through many altitude levels each with its own wind speed and direction, resulting in the satellite observation of a distorted plume that may not be well-represented by an single local wind vector. These issues further complicate the use of wind data in reconstructing the time-averaged behavior.

Here we have made novel modifications in adopting this general approach which do not rely on prior knowledge of the winds and therefore do not inherit errors associated with NWP winds. In reconstructing the time-averaged plume, we exploit the geometry of the emitted ${\mathrm{SO}}_2$ clouds themselves to inform the wind-reconstruction. To perform the plume rotation, we modify a rudimentary computer vision algorithm designed for object detection in medical imaging to detect plume-like objects in gridded ${\mathrm{SO}}_2$ data. In light of the observation that many emissions of ${\mathrm{SO}}_2$ are offset from their sources along a different vector from their primary dispersal, we also allow short translations of the data so that these puffs can still be reconstructed in the case that there is a significant lag between emission and satellite overpass time.

In the fitting stage, we employ a solution to the general time-averaged dispersion of ${\mathrm{SO}}_2$ from a point source. This fitting function incorporates the time-averaged physics of ${\mathrm{SO}}_2$ decay, advection and turbulent mixing from a point source. This fitting function is more flexible than that employed by previous studies since it more closely approximates stream-wise eddy-diffusion mixing and is better suited to examining persistent point source emissions.

We demonstrate these techniques using ${\mathrm{SO}}_2$ data generated by a newly developed probabilistic layer height and VCD algorithm designed for the Cross-track Infrared Sounder (CrIS), a hyperspectral infrared (IR) sensor aboard the Joint Polar Satellite System (JPSS) series satellites: Suomi-NPP (SNPP) and NOAA-20 [6]. Since this algorithm retrieves the ${\mathrm{SO}}_2$ layer height, the retrieved VCD estimate is already calibrated for the appropriate altitude. This CrIS ${\mathrm{SO}}_2$ algorithm is currently integrated into the NOAA/CIMSS VOLcanic Cloud Analysis Toolkit VOLCAT; [7,8,9,10], where it is used to generate ${\mathrm{SO}}_2$ alerts as well as probabilistic cloud object properties for real-time detection, characterization, and tracking of volcanic clouds in support of aviation safety.

Although the use of hyperspectral IR ${\mathrm{SO}}_2$ data limits the accuracy of plume measurements in the lower troposphere, particularly in very moist tropical atmospheres, this data can still make useful flux estimates at high-latitude volcanoes and at low-latitude, but high-altitude volcanoes. Of particular importance, IR ${\mathrm{SO}}_2$ data can fill an important data gap in the UV-based record: estimating ${\mathrm{SO}}_2$ emissions from high-latitude volcanoes through the polar winters when there is insufficient solar backscatter for UV sensors due to the high solar zenith angle (SZA).

Our approach is modular, containing three main components: (1) the gridded ${\mathrm{SO}}_2$ data including a height estimate, (2) a wind-direction reconstruction algorithm to use in constructing time-averaged plumes, (3) an independent estimation of one of the sought geophysical parameters (here, ${\mathrm{SO}}_2$ decay rate or lifetime) and (4) a fitting model to estimate the plume shape and source flux. Consequently, these components are interchangeable with only minor adjustments which allows for customized application of this technique to suit the circumstances of measurement. Although we demonstrate this technique here with IR data, it is sufficiently general to be applied to UV measurements provided that an ${\mathrm{SO}}_2$ layer height can be reliably estimated. Additionally, these methods could be incorporated into a hybrid UV-IR technique such as changing from UV- to IR-based ${\mathrm{SO}}_2$ in the winter at high-latitude volcanoes.

1.1. Basic Principles of Time-Averaged SO${}_2$ Flux Estimation

As described above, similar techniques for estimating time-averaged ${\mathrm{SO}}_2$ flux ($\dot{m}$) rely on the ability to reconstruct a representative plume for a given period of time, requiring that the atmospheric conditions, especially wind direction, are remarkably consistent, or that separate analyses can be partitioned by wind direction and speed, or that each plume observation can be rotated as if all were advected by a constant wind. Once a representative time-averaged plume is reconstructed, the flux can be fit to the VCD or in one-dimension as a line density. These fitting functions are meant to represent the time-averaged dispersion of ${\mathrm{SO}}_2$ from a point source and generally reflect some kind of spatial decay due to advection combined with turbulent diffusive mixing and photochemical oxidation kinetics. In general, the time-averaged flux is calculated as the ratio of the time-averaged cloud mass and the ${\mathrm{SO}}_2$ decay timescale or lifetime e.g., [1,2]:

$$\dot{m}=\frac{M}{\tau }.$$

(1)

Since satellite observations have direct access to the cloud mass (M), the ${\mathrm{SO}}_2$ lifetime ($\tau$) is the principal uncertain parameter and the main target of fitting. In most studies, this lifetime is interpreted in the context of first-order decay kinetics where the rate of decay of the ${\mathrm{SO}}_2$ mass in a discrete puff is proportional to the puff mass ($dM/dt=-kM$) with the apparent decay rate $k=1/\tau$. Since the time-averaged plume represents a steady state condition, it is useful to convert temporal decay into spatial decay using the mean wind speed (u) and fit the decay of the time-averaged ${\mathrm{SO}}_2$ cloud over the length scale $u/k$, using a wind dataset to estimate u.

In studies without wind-rotation, e.g., [1,11,12], the fitting function is an Exponentially Modified Gaussian (EMG) which combines downwind decay with wind-direction spreading as the convolution of an exponential and a Gaussian distribution for the ${\mathrm{SO}}_2$:

$$\mathrm{EMG}\left(x\right)=\left(\mathrm{Exp}(k/u)\ast \mathcal{N}(0,{\sigma }^2)\right)\left(x\right)$$

(2)

where $\sigma$ is an empirical downwind spreading length. Although these studies did not use wind-rotation, they used wind datasets to constrain the mean wind, and thus estimate the ${\mathrm{SO}}_2$ lifetime ($\tau =1/k$) by fitting the spatial decay of the time-averaged ${\mathrm{SO}}_2$ cloud.

In studies using wind-rotation, the fitting model proposed by [2] and implemented by [3,4,5,13] is the product of a traditional one-dimensional EMG model along the wind direction and a dispersion solution from classical Gaussian plume theory which allows crosswind spreading only. This approach enables the fitting function to spread out laterally downstream in a physically realistic way and still accommodate some mass upstream of the source due to the EMG. Ref. [2] introduced an empirical relationship between the two spreading parameters (related to the downwind and crosswind eddy diffusivities); however, this fitting function still requires the fitting of the wind speed, the lifetime ($1/k$), the flux, and at least one measure of spread (related to eddy diffusivity). The fitting function we introduce here is no different in regards to the number of parameters; however, it is a direct solution to the equation governing the theoretical dispersion of ${\mathrm{SO}}_2$ detailed below. In both the previous work and this contribution, the set of fitting parameters do not form an independent basis and thus one must be estimated independently. In these previous works, ancillary wind data was used to perform the reconstruction and thus a mean wind speed estimate (u) was readily available. Here, we use no ancillary datasets and choose instead to estimate the ${\mathrm{SO}}_2$ decay rate (k) within the averaging window as detailed below.

1.2. CrIS SO${}_2$ Data

In this study, we demonstrate our analysis technique using ${\mathrm{SO}}_2$ data from the hyperspectral IR CrIS instruments aboard SNPP and NOAA-20 satellites. CrIS is a Fourier transform Michelson interferometer, having a local time ascending node (LTAN) of 1:30 PM with NOAA-20 operating approximately 50 min ahead of SNPP [14]. The orbital offset of SNPP and NOAA-20 allows equatorial swath gaps in ${\mathrm{SO}}_2$ data from each CrIS instrument to filled by data from the other. The SNPP and NOAA-20 CrIS instruments have generated infrared spectra that include the ${\nu }_3$ (1300–1410 cm${}^{-1}$) ${\mathrm{SO}}_2$ absorption band at a spectral resolution of 0.625 cm${}^{-1}$ since December 2014 for SNPP CrIS (except 26 March–1 August 2019) and since February 2018 for NOAA-20 CrIS. CrIS scans between $\pm {48.3}^{°}$ with each scan containing 30 fields of regard (FOR), each of which contains 9 circular fields of view (FOV) arranged as a $3\times 3$ array. As the scanning mirror moves, the FOR rotates slightly and consequently, FOVs transition from 14 km-diameter circles at nadir, to partially overlapping 43.6 km × 23.2 km (major and minor axes) ellipses on the edge of scan leading to a swath width of approximately 2200 km with approximately 30% gaps and an average sampling distance of 16 km [14,15]. For this analysis, we converted the un-gridded FOV-level ${\mathrm{SO}}_2$ data into 16 km-pixel grids by nearest neighbor interpolation.

The data itself is generated by a newly developed probabilistic layer height and VCD algorithm which retrieves probability distributions for ${\mathrm{SO}}_2$ layer height and VCD from IR spectra in the ${\nu }_3$${\mathrm{SO}}_2$ absorption band [6]. The CrIS ${\mathrm{SO}}_2$ algorithm is a probabilistic modification of the trace-gas method of [16] to detect ${\mathrm{SO}}_2$ layer altitude and VCD. The novel aspect of this method is that uncertainty related to spectra with trace amounts of (or without) ${\mathrm{SO}}_2$ is propagated into uncertainty about the ${\mathrm{SO}}_2$ layer height and VCD, each of which is given by a probability density function (PDF; [6]). The PDF for the layer height informs the estimation of total VCD. Because this method exploits spectral residuals in the ${\nu }_3$${\mathrm{SO}}_2$ absorption band, it is largely insensitive to the presence of volcanic ash and hydrometeors as long as the ${\mathrm{SO}}_2$ is not completely obstructed [17]. Although there is considerable interference from water vapor in this band, it is exactly this interference that allows for robust height estimation owing to the significant variation of water vapor with altitude and the sharp absorption lines imparted to the spectrum by the presence of ${\mathrm{SO}}_2$ [6,16]. Despite the theoretical limitation posed by the degradation of sensitivity at very low altitudes where there is significant water vapor (e.g., in the tropical lower troposphere), this technique has thus far successfully detected and characterized many low-altitude ${\mathrm{SO}}_2$ clouds. By design, this algorithm can capture the ${\mathrm{SO}}_2$ present in a wide variety of background atmospheres. Because CrIS does not rely on backscattered sunlight for measuring ${\mathrm{SO}}_2$ and is sensitive to ${\mathrm{SO}}_2$ in a wide variety of background atmospheres, we use all available CrIS measurements in this effort without regard to SZA or cloud fraction as must be considered when using UV-based data. The only data omitted are measurements from SNPP CrIS FOV 7 which has unfavorable noise characteristics [6]. Because of these factors, each CrIS instrument provides nearly global coverage twice per day. In this study, we combine data from both CrIS instruments to increase the frequency of measurements, increasing the signal to noise ratio of the time-averaged plume reconstructions. At high-latitudes where there is significant swath overlap between orbits, we use data from multiple orbits to generate a composite snapshot for a region surrounding the target volcano. Using all data from both CrIS instruments yields at least twice daily coverage for equatorial volcanoes and approximately 4-times daily coverage for most other volcanoes.

In this work, we use a statistical measure of ${\mathrm{SO}}_2$ column enhancement (z-score), which is native to the CrIS retrieval, to define the plume geometry and wind reconstruction. Throughout this work, except where noted, we abbreviate the z-scores as z. In principle, if another source of ${\mathrm{SO}}_2$ data were used (e.g., UV spectrometers such as the Ozone Monitoring Instrument (OMI), the Ozone Mapping and Profiler Suite (OMPS), or the TROPOspheric Monitoring Instrument (TROPOMI)) the ${\mathrm{SO}}_2$ VCD values themselves could be used directly along with a critical threshold indicating significant enhancement above the background. Although we are interested in persistent ${\mathrm{SO}}_2$ sources with a typically weak signal, the CrIS-based retrieval of ${\mathrm{SO}}_2$ VCD as in [6] is most reliable for individual measurements with $z\ge 5$. However, reasonable, but noisier VCD retrievals have been accomplished down to a threshold $z\ge 3$. Although physically meaningful spatial patterns of z-score often emerge even when $z<3$, the layer height and VCD are only retrieved here for $z>3$ to ensure that the retrieval algorithm is computationally efficient and less prone to noisy results. In this application, we use the more reliable measurements to fill in areas of the swath where the VCD was not retrieved directly. Because of the averaging used in this type of analysis, the errors that this generates largely cancel out. For each gridded overpass of the target volcano’s surroundings, we compute a conversion factor, referred to here as the “gain” (g), which converts between z-scores and VCD obtained from the full probabilistic retrieval, $g=X/z$, for all measurements with $z>3$. We then estimate the geometric mean gain value ($\widehat{g}$) from these more reliable samples (those with $z>3$). We have chosen the geometric mean here since it is a more robust estimator for ratios. If no samples with $z>3$ are available, the overpass is assumed to contain only low-amplitude noise and the mean gain and VCD are set to zero. We use the mean gain value to construct an estimated VCD grid as $\widehat{X}=\widehat{g}z$. By contrast to simply using the original VCD values, the new VCD estimate has the same distribution shape as the z-scores and when multiple grids are combined, better reduces background noise than using the raw VCD which is only estimated for samples with $z>3$. This approach is well-suited to analysis of long-term ${\mathrm{SO}}_2$ flux because individual measurements are of little interest and instead, the aggregate behavior is sought. Additionally, the gain-based VCD estimation incorporates altitude information without estimating a single mean altitude, instead constructing a representative scaling between the z-scores, which contain no altitude information, and VCD values, which do contain altitude information. This approach allows the plume information below $z=3$, which can be significant in a weak plume, to be used in reconstructing the time-averaged plume without introducing large errors from less reliable VCD retrievals.

2. Materials and Methods

2.1. Object Detection-Based Mean Plume Construction

To generate a time-averaged plume that can be used to compute time-averaged geophysical quantities, including the source flux, sets of observations of a plume must be aggregated. First, raw ${\mathrm{SO}}_2$ data must be aggregated to form a representative “snapshot" of the real plumes on an azimuthal projection grid. This requires stacking up a sufficient number of overpasses to span the study domain in question. For a small domain around a given volcano, a sufficient snapshot can be generated in as little as one overpass; however, for large domains, several precessing orbits might be necessary to span the domain of study. The number of orbits depends upon the domain size, the satellite’s mean motion, and the latitude of the volcano. Once a time series of snapshots has been constructed, they must be further aggregated to construct a time-averaged empirical plume model. This empirical model represents the time-averaged steady-state plume if the wind were blowing in a constant direction, thus requiring rotation and regridding of each plume snapshot before stacking the data. Additionally the longer the time interval over which the average is constructed, the more plume snapshots are included and consequently the signal-to-noise ratio increases. In particular, for the probabilistic CrIS ${\mathrm{SO}}_2$ retrieval, ${\mathrm{SO}}_2$ total column mass loading is a normal random variable and the central limit theorem predicts that the signal-to-noise ratio (ratio of mean and standard deviation) grows as the square root of the number of repeat images. Consequently, the use of both SNPP and NOAA-20 CrIS sensors yields an approximate $40\%$ increase in signal to noise ratio over the use of only one of these two sensors.

2.1.1. Plume Detection, Source Reconstruction, and Wind Rotation

Recent efforts to characterize volcanic and industrial contaminant sources have utilized NWP or reanalysis wind data sets to rotate plume data into a common wind direction for improved analysis of persistent source fluxes [2,3,4,5,18]. Here we present an approach that does not use ancillary data sets, instead, using only image analysis of the snapshots themselves. Our technique is based on identifying near-source “plume-like” ${\mathrm{SO}}_2$ cloud objects in these snapshots and then rotating and translating them such that they originate from the volcano and are elongated in a common wind direction. This effort is considerably simpler than traditional object detection for computer vision since the input data are gridded ${\mathrm{SO}}_2$ measurements and therefore lack other types of image features that might interfere in a more general image processing application. In particular this method estimates a dominating wind direction by searching the space of all lines through the ${\mathrm{SO}}_2$ z-score image from all azimuthal projections to find the projection angle giving the largest total line integral. This corresponds roughly to finding the azimuthal angle which looks down the elongation axis of a plume. Once this projection angle is found, we generate a normalized profile of the data along the plume in the projection direction and estimate the approximate location of the suspected source along this profile. This is generally not the true source of the ${\mathrm{SO}}_2$, but instead attempts to approximate the location to which the source region of the plume has been advected. Because this special projection direction is taken to be the wind direction, the grid is rotated so that the wind is in a specified direction, and the grid is translated so that the estimated source region is moved to the origin. In this scheme, all emissions observations undergo the same processing including more equant “blob-like” emissions that do not conform to simplified plume theories. These emissions are assumed to be the result either of an instantaneous puff of ${\mathrm{SO}}_2$ advecting downwind or more continuous emissions which spread nearly equally in all directions, indicative of the weak role of downwind dispersal compared with turbulent mixing. Because even a short duration, but continuous emission will tend to stretch out downwind if the mean wind is strong relative to turbulent mixing, in this work, we implicitly assume that all observed emissions are of this latter type.

Mathematically this is accomplished by analyzing the Radon transform, used widely in medical tomography, signal processing of seismic data, and computer vision e.g., [19,20,21,22]. In this formulation, although it is ${\mathrm{SO}}_2$ vertical column density (VCD) that will be used to construct the time-averaged plume, we use the CrIS ${\mathrm{SO}}_2$ z-scores (Figure 1a) to define the necessary grid transformations owing to the beneficial statistical properties of the z-score compared with the CrIS-based VCD retrieval for weak ${\mathrm{SO}}_2$ plumes [6,16,23]. If UV-based ${\mathrm{SO}}_2$ data were used, the VCD estimate itself could be used instead.

For a gridded snapshot of ${\mathrm{SO}}_2$ z-score ($z(x,y)$), we define the Radon transform on the set of all lines in the plane as

$$\mathcal{R}z(\alpha ,s)={\int }_{\mathcal{L}(\alpha ,s)}z(x,y)\mathrm{d}t$$

(3)

where each projection line $\mathcal{L}(\alpha ,s)$ has been parameterized as

$$\left\{\begin{array}{ll}x(t;\alpha ,s)=t\sin \alpha +s\cos \alpha & \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{(4a)}\\ y(t;\alpha ,s)=-t\cos \alpha +s\sin \alpha & \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{(4b)}\end{array}\right.$$

with t, an arc-length parameter, measuring true distance along the line. Here, s is the perpendicular distance from the origin to the projection line and the angle $\alpha$ is the angle from the positive x-axis to the projection line normal vector $(cos\alpha ,sin\alpha )$ which is also the angle between the positive x-axis and the line segment of length s perpendicular to $\mathcal{L}$ (Figure 1c). The tangent vector of $\mathcal{L}(\alpha ,s)$ is $(sin\alpha ,-cos\alpha )$. From the transformed map $\mathcal{R}z(\alpha ,s)$ or “sinogram" (Figure 1d), we define the special point

$$(\widehat{\alpha },\widehat{s})=\underset{(\alpha ,s)}{arg\; max}\mathcal{R}z$$

(5)

which corresponds to the line $\widehat{\mathcal{L}}=\mathcal{L}(\widehat{\alpha },\widehat{s})$ that most reflects an elongated plume in the image with the projection vector $\widehat{v}=(sin\widehat{\alpha },-cos\widehat{\alpha })$. We then construct a profile of the mapped z-scores along this line:

$$\widehat{z}\left(t\right)=z(x(t;\widehat{\alpha },\widehat{s}),y(t;\widehat{\alpha },\widehat{s})).$$

(6)

Because this profile is taken along a straight line, it does not accurately capture the plume axial profile since the plume axis may be curved or distorted by the flow as in the subtle von Kármán vortex street distorting the plume axis in Figure 1 and Figure 2. To account for this, we search the space of nearby profile lines that are parallel to $\widehat{\mathcal{L}}$ and extract the maximum value of the z-score for each value of t across all of these lines. This gives an approximation to the plume axial profile, though it is slightly foreshortened locally where the real plume axis makes excursions from the line $\widehat{\mathcal{L}}$. Once the overall profile is constructed, negative values are removed and the profile is normalized (Figure 2a). We then compute the 10th, 50th, and 90th percentiles of this modified profile, respectively. Based on the asymmetry of the 10th and 90th percentiles about the median, we accept or reverse the projection direction to determine the plume dispersal direction, reflecting the assumption that the plume is most concentrated near the source and is more dispersed in the downwind direction, that is, the profile is skewed downwind (Figure 2b).

With the projection direction now matched to the plume dispersal direction, the plume source must be estimated. As mentioned above, the approximate source region is meant to describe the location of the source region relative to the plume’s current configuration, meaning it may have been advected some distance away from the volcano after emission by the time of imaging. To estimate this position along the profile, $t_{source}$, we exploit the simple one-dimensional advection-diffusion-reaction model (detailed below) which we use to estimate $t_{source}$ based on an estimate of the Peclet number and upstream-downstream mass partitioning for the profile (Figure 2c,d). This technique is detailed in Appendix A.

From these estimates, the source position is taken as $(t_{source},\widehat{\alpha },\widehat{s})$ in the projection coordinates and the Cartesian coordinates of the source position can be found from Equation (4a,b). The Cartesian grid coordinates are then translated so that the source is at the origin and rotated so that the dispersal direction is in the pre-specified direction (Figure 2e,f).

2.1.2. Filtering and Smoothing in Wind Reconstruction

Within this plume rotation and source reconstruction technique, several additional steps are taken to improve the quality of the reconstruction. At multiple stages detailed below, we make use of exponential smoothing kernels (Appendix B) rather than Gaussian smoothing or other types of smoothing for two reasons: (i) exponential smoothing reduces noise variance to a greater extent than Gaussian smoothing [24], and (ii) the exponential dependence in our chosen physics-based model suggests that exponential smoothing will only minimally distort the edges of the ${\mathrm{SO}}_2$ plume data.

First an attempt is made to remove background ${\mathrm{SO}}_2$ which may come from a bias in the noise or from anthropogenic ${\mathrm{SO}}_2$ sources. Since the volcanic ${\mathrm{SO}}_2$ plume occupies such a small fraction of a given grid, we first subtract the grid mean value from all pixels. Next, we filter the grid with a special filter designed to heavily weight measurements that are close to the volcano and are also in a neighborhood of the strongest detections (Figure 1b). This two-part filter is constructed as follows.

First, a radial weighting filter is constructed from a truncated radial Gaussian function centered on the volcano in question:

$$\mathcal{G}\left(r\right)=\left\{\begin{array}{cc}exp\left(-r^2/2{\sigma }_r^2\right)\hfill & r<r_{max}\hfill \\ 0\hfill & r\ge r_{max}.\hfill \end{array}\right.$$

(7)

This weighting filter has unit weight at the volcano and is zero outside of the radius $r_{max}=4{\sigma }_r$. This multiplicative filter is designed specifically to heavily weight plumes within ${\sigma }_r$ (the standard deviation of the Gaussian) from the volcano. This filters out ${\mathrm{SO}}_2$ that has drifted into the edges of the image from other sources as well as plumes from the target volcano that were emitted previously to limit the possibility of double counting. We use this radially decaying weighting filter instead of a fixed search radius because it allows for variations in the wind speed which may advect puffs of ${\mathrm{SO}}_2$ to variable distances away from the volcano.

The second filter component generates a buffer around all z-score values above a threshold value ($z_{th}$). This is accomplished by generating a binary map of all z-scores above the threshold and then convolving the binary map with an exponential kernel (K) with 20-km e-folding radius. A binary image is then generated to indicate where the numerical convolution result is nonzero and then this new binary image is again convolved with the same exponential kernel. The resulting map is of unit weight over a buffer patch very close to the high-confidence detections and decays to zero over a distance of approximately 50–100 km. This buffering weight ($\mathcal{B}$) can be summarized mathematically as

$$\left\{\begin{array}{ll}\mathcal{B}=K\ast {\mathbf{1}}_{\mathsf{\Omega }}& \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{(8a)}\\ \mathsf{\Omega }=\{K\ast {\mathbf{1}}_{\{z\ge z_{th}\}}>0\}\hfill & \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{(8b)}\end{array}\right.$$

where ${\mathbf{1}}_{\{·\}}$ is the indicator function for the specified region and the convolution is understood to be performed numerically with a compactly supported kernel K. In our modular framework, if UV-based ${\mathrm{SO}}_2$ data were used instead of the IR data used here, the VCD estimate itself along with a threshold value for significance would be used to construct this filter.

An overall weight ($\mathcal{W}$) is constructed as the point-wise product of the detection buffer weight and the radial weighting filter ($\mathcal{W}(x,y)=\mathcal{G}(x,y)\mathcal{B}(x,y)$, Figure 1b). As a last step, we multiply this overall weighting map $\mathcal{W}$ by the original z-scores and perform a minimal smoothing by convolution with another exponential kernel of only 4-km e-folding radius to generate a smoothed and filtered z-score image which can be more readily analyzed by the above method to find plume-like ${\mathrm{SO}}_2$ clouds near the volcano (Figure 1c).

This filtered and smoothed z-score image is only used to define the necessary grid transformations to reconstruct the plume as detailed above. Once this is complete, the raw z-score image is interpolated onto the new grid (Figure 2g). A final filtering and smoothing step is performed after rotation and re-gridding. This step involves the binary image segmentation (specifically, connected-component labeling) of the buffered detection region $\mathsf{\Omega }$ computed for the rotated and re-gridded z-scores. Although we use the Python-based SciPy multidimensional image processing package, similar algorithms are available in most modern code libraries. From this, we keep only the connected detection region which contains the re-gridded origin (plume source). All other pixels are set to zero and the result is smoothed with a 4-km e-folding radius exponential kernel (Figure 2g).

2.1.3. Time-Averaging and Post-Processing

After grid transformation as described above, we construct the estimated VCD for each plume snapshot and average together all available snapshots in the time-averaging window. As a final step, the time averaged grid is translated so that the maximum time-averaged VCD value is at the volcano. This process is meant to reconstruct an empirical model of the plume from the volcano under conditions of constant emission, wind field, mixing, oxidation, and dilution.

2.2. Physics-Based SO${}_2$ Cloud Fitting

Here we present a very simple physical model for the continuous point release, dispersion, and evolution of a cloud of ${\mathrm{SO}}_2$ subject to a constant wind field which we will use to throughout this work to estimate the volcanic ${\mathrm{SO}}_2$ flux. The proposed model and its solutions are well known in the literature and can be found in or synthesized from the elements of many texts on environmental fluid mechanics e.g., [25,26,27].

We model the continuous point release of ${\mathrm{SO}}_2$ and incorporate constant advection, Fickian turbulent diffusion, and first-order reaction kinetics representing total loss of ${\mathrm{SO}}_2$ from the cloud as an advection-diffusion-reaction Equation (ADRE):

$$\frac{\partial C}{\partial t}+u\frac{\partial C}{\partial x}=D_x\left(\frac{{\partial }^{2}C}{\partial x^2}+\frac{1}{{\alpha }_y}\frac{{\partial }^{2}C}{\partial y^2}+\frac{1}{{\alpha }_z}\frac{{\partial }^{2}C}{\partial z^2}\right)-kC+S$$

(9)

where u is the (constant) wind speed, k is a first-order reaction-kinetic decay rate, $D_x$ is the downstream eddy diffusivity, ${\alpha }_y=D_x/D_y$ and ${\alpha }_z=D_x/D_z$ are the aspect ratios comparing the downstream eddy diffusivities to those in the crosswind horizontal and vertical directions and S is a source term describing continuous point release of ${\mathrm{SO}}_2$ from the origin with mass flow rate $\dot{m}$. In this section, z refers to a vertical coordinate and not ${\mathrm{SO}}_2$ z-scores used elsewhere in the text. Note that as in [28] the chemistry (decay) term $-kC$ represents the ${\mathrm{SO}}_2$ loss due to actual processes (oxidation and conversion to sulfates, dry and wet deposition) as well as apparent ${\mathrm{SO}}_2$ loss related to dilution of the cloud below the sensor detection threshold. In this very simplified model, we have assumed that all eddy (turbulent) diffusivities are Fickian, that is, they are constants for a given analysis, or alternatively, they change sufficiently slowly that differences are not relevant over the time period of analysis. Other than the diffusivity aspect ratios which directly control the cloud’s spread, the shape of the cloud is also governed by the Peclet number of the flow which is typically defined as the ratio of advective transport to diffusive transport. Since the problem is set in an unbounded domain (${\mathrm{SO}}_2$ dispersion from an elevated point source in free atmosphere), there is no single characteristic length. In this study we choose as a characteristic length, the “e-folding distance” ($u/k$)—the distance traversed in one e-folding timescale ($\tau =1/k$) under constant advection. This length scale, the distance over which a puff of ${\mathrm{SO}}_2$ decays by a factor of e in constant wind, is the primary measurement that previous studies have used to estimate the ${\mathrm{SO}}_2$ lifetime given knowledge of wind speed e.g., [1,2,4,18]. We therefore define the Peclet number for the problem as

$$\mathrm{Pe}=\frac{u^2}{D_{x}k}.$$

(10)

This Peclet number controls the ${\mathrm{SO}}_2$ plume’s elongation in the wind direction and the upstream-downstream distribution of mass.

2.2.1. Application of Simplified Theory to Satellite Measurements of Volcanic ${\mathrm{SO}}_2$  Plumes

Although the governing ADRE can be used to describe a wide variety of elementary flows, we focus on analysis of the cloud produced by a continuous point release of ${\mathrm{SO}}_2$ in a time-averaged sense. As long as the time-averaging interval is much larger than characteristic changes in the ${\mathrm{SO}}_2$ cloud between the time interval endpoints, the time-averaged ADRE is asymptotically equivalent to the steady state ADRE with a constant source flux. We focus on this steady state solution in the present work since we will be analyzing time-averaged ${\mathrm{SO}}_2$ plume data integrated over time intervals much longer that the timescales of meaningful cloud dynamics predicted by the transient model (Appendix C). Additionally, from the perspective of a satellite, the above three-dimensional model is not very helpful for a trace gas like ${\mathrm{SO}}_2$ for which there is little or no vertical profile information. Consequently, most space-based ${\mathrm{SO}}_2$ measurements are of VCD, denoted here as $X(x,y)$. This measurement destroys any dependence on the vertical diffusivity (and consequently the parameter ${\alpha }_z$) since the VCD is a quantity integrated over the entire vertical coordinate space after [27]:

$$X(x,y)=\int C(x,y,z)\mathrm{d}z=\frac{\dot{m}\sqrt{{\alpha }_y}}{2\pi D_x}exp\left(\frac{ux}{2D_x}\right)K_0\left(\frac{v\rho }{2D_x}\right)$$

(11)

where the following shorthand has been used:

$$\left\{\begin{array}{ll}v=\sqrt{u^2+4D_{x}k}& \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{(12a)}\\ \rho =\sqrt{x^2+{\alpha }_y{y}^2}& \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{(12b)}\end{array}\right.$$

and $K_0$ is the zero-th order modified Bessel function of the second kind. Here v is the average velocity of transient perturbations (Appendix C) and describes the excess dissipation caused by reaction and turbulent diffusion processes and $\rho$ is an elliptical radius measuring contours of constant dispersion in the absence of wind. Of significant importance for model fitting, the parameters in this equation do not form an independent basis. In particular, there is not unique dependence on the downstream eddy diffusivity ($D_x$). To clarify this point, we rewrite this solution in terms of independent parameters only:

$$X(x,y)=\frac{\mu \sqrt{{\alpha }_y}}{2\pi }exp\left(\frac{\omega x}{2}\right)K_0\left(\frac{\lambda \rho }{2}\right)$$

(13)

where $\mu =\dot{m}/D_x$, $\omega =u/D_x$, and $\lambda =v/D_x=\sqrt{{\omega }^2+4\kappa }$ with $\kappa =k/D_x$.

This VCD model, equivalent to the solution for an infinite column source at the origin, depends only on the parameters ${\alpha }_y$, $\mu$, $\omega$, and $\kappa$. If space is re-scaled by the characteristic length (i.e., ${\mathbf{x}}^{*}=\mathbf{x}k/u$), it is easily demonstrated that shape of this cloud is governed only by the crosswind aspect ratio and the Peclet number since $\omega x=\mathrm{Pe}{x}^{*}$ and $\lambda \rho =\sqrt{{\mathrm{Pe}}^2+4\mathrm{Pe}}{\rho }^{*}$.

2.2.2. Relationship to Previous Fitting Model

By contrast to the ADRE solution used here, the fitting function of [2] can be written as $X(x,y)\propto \tilde{X}(x,y)·\mathrm{EMG}\left(x\right)$ where the Gaussian plume solution ($\tilde{X}$) approximately solves the steady ADRE in the high Peclet limit in which stream-wise diffusion and decay are neglected:

$$u\frac{\partial \tilde{X}}{\partial x}=D_y\frac{{\partial }^2\tilde{X}}{\partial y^2}+\dot{m}{\delta }^2(x,y)$$

(14)

As described earlier, this approach from [2] enables the fitting function to accommodate some mass upstream of the source and to spread out downstream although it is heuristic, lacking a direct solution from the physics of the ADRE. The 2D EMG model has the property that its peak is downwind of the source. This could be a useful property if the average lag time between emitted pulses and satellite overpass could be estimated. Since this is not known, our plume source-reconstruction approach is more flexible since it more closely approximates upstream turbulent mixing and is better suited to examining persistent point source emissions. Additionally, this fitting model has exactly the same issues related to fitting non-uniqueness as the ADRE used here, requiring independent estimation of one of the geophysical parameters, that is, the flux, wind speed, decay rate, and measures of spread are sought, but one of these must be estimated independently. In previous studies using wind data, the wind speed is estimated and the other parameters are fit.

2.2.3. Fitting the ADRE Parameters

Although the VCD model $X(x,y;{\alpha }_y,\mu ,\omega ,\kappa )$ could itself be fit to time-averaged satellite ${\mathrm{SO}}_2$ VCD observations in a similar manner as in [2,4,5,12], the signal to noise ratio can be improved by another integration, specifically integrating through the crosswind (y) direction as in [1,11,12,18] to generate a line density. Although this will destroy information related to the crosswind aspect ratio ${\alpha }_y$, this has no effect on measuring the distribution of ${\mathrm{SO}}_2$ mass upstream and downstream which is of greater consequence than measuring the aspect ratio. Additionally, once the line density parameters are fit, the aspect ratio can be estimated by other means (Appendix D). Because the main parameters of interest to the present work are estimates of the flux $\dot{m}$, the ${\mathrm{SO}}_2$ lifetime $1/k$, and to a lesser extent the Peclet number of the flow, integrating out the crosswind dependence will not significantly impact the estimation of these metrics. In practice, we integrate over the $\pm 2$ standard deviation width of the time-averaged plume which we determine by first integrating in the downwind direction and then treating as a probability distribution. Performing this integration gives the profile of ${\mathrm{SO}}_2$ mass per distance downwind contained in the plane perpendicular to the plume axial line or wind vector:

$$\mathcal{L}\left(x\right):=\int X(x,y)\mathrm{d}y=\frac{\mu }{\lambda }exp\left(\frac{\omega x-\lambda \left|x\right|}{2}\right).$$

(15)

As VCD refers to the density of the species contained in a vertical column, we refer to this quantity (a line density of ${\mathrm{SO}}_2$ in planes perpendicular to the plume axis) as the axial planar density (APD).

The one-dimensional APD solution is equivalent to the solution to the three-dimensional governing equation for an infinite planar source at $x=0$. This integration provides four significant benefits: (i) APD data have a stronger signal to noise ratio than VCD data yet preserves the same parametric dependence on $\mu$, $\omega$, and $\kappa$ as $X(x,y)$, (ii) The use of APD instead of VCD eliminates the possibility of model error due to ${\mathrm{SO}}_2$ data asymmetries in the crosswind direction, (iii) Regardless of the crosswind eddy diffusivity aspect ratio, the upstream-downstream ${\mathrm{SO}}_2$ mass partitioning will not be affected so an accurate Peclet number can still be computed from APD fitting and (iv) the high signal-to-noise ratio and simpler functional form make nonlinear fitting of $\mathcal{L}\left(x\right)$ to APD data more likely to converge owing to the ease of generating high-fidelity initial parameter estimates compared with what is required for fitting $X(x,y)$ to VCD data.

Because of these benefits, we construct an APD profile from the time-integrated plume VCD and then fit to it the model $\mathcal{L}(x;\mu ,\omega ,\kappa )$. With these fitting parameters, we can immediately calculate the plume Peclet number as $\mathrm{Pe}={\omega }^2/\kappa$. An important benefit of this model is that functionally, it is an exponential decay upstream and downstream with a matched point at the source. Exploiting this simplicity, we fit a log-transformed APD profile rather than the APD profile itself since this is closely related to linear regression up- and downstream with approximately Gaussian error, enabling the use of traditional fitting metrics such as the coefficient of determination ($R^2$) without incurring the significant difficulties in interpreting nonlinear regressions this way e.g., [29].

To estimate the degassing flux ($\dot{m}$), we must independently estimate one of u, k, $D_x$, or $D_y$. If, as in similar studies e.g., [1,2,3,4,5,11,12,18], wind data is already used to construct the time-averaged plume, then the mean wind speed could be used compute $(\dot{m},D_x,k)$ from $(\mu ,\omega ,\kappa )$. Here, we focus on the decay rate k since to first order, it can be estimated directly (with uncertainty) from the record of individual plume snapshots that make up the time-averaged composite. Once this is estimated, the quantities of interest ($\mathbf{q}$) are determined from the parameter estimates ($\mathbf{p}$) by

$$\mathbf{q}=\frac{k}{\kappa }\mathbf{p}$$

(16)

where the quantities of interest are $\mathbf{q}={(\dot{m},u,D_x,D_y)}^{⊺}$ and $\mathbf{p}={(\mu ,\omega ,1,{\alpha }_y^{-1})}^{⊺}$.

2.2.4. Decay Rate Estimation

From [6], the CrIS ${\mathrm{SO}}_2$ total mass (M) is a normal random variable as is the time derivative $dM/dt$. Consequently, if we restrict analysis to periods with $dM/dt<0$, the decay rate (k) and the ${\mathrm{SO}}_2$ lifetime ($\tau =1/k$) are both ratios of normal random variables with an almost certainly non-zero variable as the denominator. Ratio distributions of this type can be approximated as being lognormally distributed [30]. The distributions of estimated ${\mathrm{SO}}_2$ lifetimes from major sources globally in [3] exhibits an approximately lognormal shape, lending evidence for the implementation of a lognormal model. We use this fact to compute a lognormal PDF for likely values of the decay rate for each time-averaged plume and propagate that uncertainty to construct a similar PDF for the mass flux and other variables as well.

Since a lognormal distribution is parameterized by the Gaussian parameters of the log-transformed variable, we estimate the parameters of the lognormal distribution using the valid log-transformed decay rate samples. We generate our decay rate samples from the time series of plume observations in the time-integration interval using the simplifying assumption that whenever the cloud mass in the domain is decreasing, the flux of ${\mathrm{SO}}_2$ from the volcano at that time is much less than the magnitude of the mass loss rate. Such an assumption is related to treating the time-averaged plume as the average of a series of discrete puffs. This is a reasonable assumption for lower-middle tropospheric ${\mathrm{SO}}_2$ injection where the ${\mathrm{SO}}_2$ lifetime is generally regarded as being significantly shorter than for major eruptions injecting ${\mathrm{SO}}_2$ into the upper troposphere or stratosphere. Since smaller, less-powerful ${\mathrm{SO}}_2$ emission is the target activity, the conditions of this assumption appear reasonable.

This approach yields samples of the decay rate as

$$k_i=\frac{1}{M_i}\frac{M_i-M_{i+1}}{t_{i+1}-t_i}$$

(17)

for all i where $M_{i+1}<M_i$ and $M_i$ is the i-th total ${\mathrm{SO}}_2$ cloud mass in the sequence of plume images. Note that in this scheme, the maximum value of the samples is limited by the time interval as $k_i^{max}=1/\mathsf{\Delta }{t}_i$ which always occurs when the cloud mass decays to zero (i.e., $M_{i+1}=0$, or apparent 100% loss in one time interval). We then correct the samples by a specialized procedure detailed in Appendix E which accounts for this hard upper limit (right-censoring). The corrected samples are then used to estimate the parameters of the lognormal distribution for the decay rate.

2.3. Uncertainty Quantification

From [6], CrIS VCD data generally agrees very well with high-resolution data from TROPOMI for small to moderate columns (<50 DU) and tends to underestimate very large TROPOMI columns by about 25% consistently across a variety of scales. Because the vast majority of the ${\mathrm{SO}}_2$ data analyzed here are below this threshold and TROPOMI data uncertainties are approximately 30–50% of those of OMI [13], we estimate that the errors originating from the use of CrIS to compute total cloud mass are somewhat less than those from OMI. We estimate the uncertainty on the CrIS VCDs in each plume image as approximately 14% for these relatively dilute plumes and as much as 20% for plumes with a large fraction of columns $>50$ DU, though these are uncommon outside of significant eruptive events. The estimation of VCD from z-scores is expected to have a very small impact on total error.

In the absence of errors due to the rotation, the uncertainty from stacking multiple images should decay as the number of stacked images increases, specifically decaying as $1/\sqrt{n}$. Ref. [3] estimated errors from the NWP wind rotation technique as approximately 6% for wind speed and direction and approximately 20% for wind height. Similar to [3], varying the detected wind direction with $1\sigma$ uncertainty of ${15}^{\circ }$ yields approximately 4–7% variation in our time-integrated plume. Additionally, our source reconstruction technique yields between 17–35% variation compared with the traditional approach without source reconstruction-based wind rotation.

Additional error is incurred due to the highly unsteady nature of volcanic emissions. Compared with anthropogenic sources, volcanic outgassing is highly variable, varying by orders of magnitude on a variety of timescales e.g., [31,32,33,34,35,36]. In our simple ${\mathrm{SO}}_2$ dispersion model, the fitting solution is the result of time integration of a dense sequence of discrete puffs, which is at least somewhat accurate to nature; however, error is incurred whenever the time delay between a new discrete puff and the next satellite overpass is maximized, that is when a discrete puff occurs just after an overpass. This error is due mainly to the chemistry and turbulent diffusion effects which both dissipate the concentration and total mass of the puff by some amount when the next overpass occurs. Because at low altitudes the ${\mathrm{SO}}_2$ lifetimes are on the order of hours, significantly less ${\mathrm{SO}}_2$ could possibly be measured in the next overpass for that puff than if the puff happened just before an overpass. In general, for a re-imaging time of T and an ${\mathrm{SO}}_2$ lifetime $\tau$, the fraction ($\phi$) of puff mass that can be observed at imaging time for a puff emitted at $t_e$ in the re-imaging interval is $\phi \left(t_e\right)=exp(-(T-t_e)/\tau )$. Assuming that puffs are uniformly distributed in each imaging window, the expected value loss for a given puff is $\mathbb{E}(1-\phi \left(t_e\right))=\{1-(\tau /T)[1-exp(-T/\tau )]\}$ by the time it is imaged for the first time. For an ${\mathrm{SO}}_2$ lifetime of 6 h e.g., [3,13] and a re-imaging window of 12 h, each puff will have decayed by an average of approximately 57%. Fortunately, this theoretical error propagates only minimally to the decay rate estimation since that is already based on relative loss, so such a consistent bias factor will play little role. If there were just a single puff released in an interval this error would be inherited in the accumulation of the time-averaged plume; however, the more puffs released, the more total mass is accumulated. In the limit of an infinite number of puffs in a re-imaging interval, the steady state fitting function would be an accurate representation of the measured distribution of ${\mathrm{SO}}_2$ and this error would be zero theoretically. Although we do not derive a theoretical correction factor to describe this effect, based on the formula above, we estimate that this loss factor is approximately $(57\pm 18)$% where the $1\sigma$ variation is derived from considering $\tau =6\pm 3$ h. This yields a uncertainty contribution from this factor of approximately 32%.

Combining these sources, we estimate that the total uncertainty on the time-averaged plume is between 37–48% over a wide range of averaging periods and only 18–36% if the effect of emission-measurement asynchrony is determined to be negligible. As described above, we fit an APD model rather than a VCD model which greatly improves the signal to noise ratio of the data. Fitting uncertainty is reduced with respect to VCD fitting and as with the time-averaging procedure, decays as inverse of the square root of the number of crosswind measurements, yielding a fitting uncertainty which is much less than these other sources.

Uncertainty in the decay rate estimation is derived mainly from uncertainty in the CrIS VCD data which is amplified somewhat by a fairly noisy decay rate estimation process (Appendix E). Even if all CrIS VCD data have the larger uncertainty of 20% and the total mass integration domain only includes about 900 measurements ($30\times 30$ pixels), the maximum uncertainty on the corrected decay rate samples is only about 10%. All remaining uncertainty on the decay rate is accounted for by representing the decay rate as a lognormal distribution. Owing to the maximum base uncertainty on the samples, we impose a minimum 10% uncertainty on the decay rate lognormal distribution.

Lastly, because the method used to estimate the decay rate is noisy, the uncertainty propagated to the quantities of interest from the decay rate (representation as a lognormal distribution with 10% minimum uncertainty) is much larger than the uncertainty from the fitting procedure and in most cases, the uncertainty from all other sources, though it is typically less than about 45% overall. Although this assumes that all errors discussed here are uncorrelated, which may be violated in some cases, the probability dispersion present in the decay rate distribution alone still likely accounts for most of the total uncertainty. For comparison, uncertainty on previous state-of-the-art annual emissions have been estimated as 55% for stronger sources (>100 kt/yr) and 67% for weaker sources (<50 kt/yr) e.g., [3,5].

Because each parameter of interest is proportional to k (Equation (16)), the retrieved PDFs of these parameters will inherit the decay rate’s lognormal distribution. Specifically, the decay rate is distributed as

$$k\sim \mathrm{Lognormal}\left({\mu }_k,{\sigma }_k^2\right)$$

(18)

where the lognormal distribution is parameterized by ${\mu }_k$ and ${\sigma }_k$, respectively, the mean and standard deviation of $lnk$. As described above, we impose a minimum uncertainty (here, coefficient of variation) for the decay rate, constraining the value of ${\sigma }_k$ as ${\sigma }_k\ge \sqrt{ln(1+{0.1}^2)}$, derived from the parameterization of the coefficient of variation for a lognormal variable. Each parameter of interest $q_i$ is therefore distributed as

$$q_i=\frac{k}{\kappa }{p}_i\sim \mathrm{Lognormal}\left({\mu }_k+ln\frac{p_i}{\kappa },{\sigma }_k^2\right)$$

(19)

and the ${\mathrm{SO}}_2$ lifetime is distributed as

$$\tau =\frac{1}{k}\sim \mathrm{Lognormal}\left(-{\mu }_k,{\sigma }_k^2\right).$$

(20)

3. Results

As described above, one of the significant benefits of developing these techniques for IR sensors is the ability to fill gaps in the UV-based record at high-latitude volcanoes in the winter when there is insufficient UV backscatter to make reliable measurements. Additionally, at high-latitude volcanoes, CrIS swaths from each JPSS satellite overlap significantly, providing at least 4 complete scans over a target volcano per day from NOAA-20 and SNPP combined. This high combined return rate and orbital coverage is ideal for increasing the signal-to-noise ratio of the time-averaged results. Below, we demonstrate these techniques for a 4-month period of unrest at Veniaminof Volcano between September–December 2018 (Figure 3).

3.1. 4-Month Summary of Eruptive Sequence

Between September and December 2018, Veniaminof volcano in the Aleutian islands underwent a period of unrest characterized by a sequence of small-moderate explosive episodes, lava flows, and a mix of passive and active ${\mathrm{SO}}_2$ emissions (Taryn Lopez pers. comm); [37]. Between September 3–4 the Alaska Volcano Observatory (AVO) raised its aviation color code (ACC) and volcano alert level (VAL) from Green/Normal to YELLOW/Advisory and then to Orange/Watch based on an increase in seismic activity and reports and webcam footage of low-level pulsatory ash emissions. On 21 November 2018, significant explosive activity sent an ash plume to a height of 15,000 ft (4.6 km), sending ash 150 km downwind which prompted AVO to raise the ACC/VAL to RED/Warning. Other than this event, AVO maintained the Orange/Watch color code and alert level through the beginning of January 2019 when the activity stopped [37]. Throughout this period, ${\mathrm{SO}}_2$ height estimates retrieved from the probabilistic CrIS algorithm were dominantly below 4 km with several plumes exceeding 5–6 km over small areas (Figure 4b). This is an ideal eruption on which to demonstrate this new flux algorithm using CrIS data since the high SZA at this latitude significantly degrades the signal-to-noise ratio of UV-based ${\mathrm{SO}}_2$ data, preventing reliable analyses of this type. Using IR data, like that from CrIS, is the only way to measure the flux consistently throughout this interval particularly between the early and late stages of this eruptive sequence.

As described above, this algorithm uses a radial Gaussian search weight parameterized by a soft search radius (${\sigma }_r$) with a hard cutoff at $4{\sigma }_r$ rather than a single hard search radius cutoff as in [3,13] and others. For the case study of Veniaminof, we used a soft search radius of 150 km, meaning that there is a hard cutoff at 600 km and that theoretically 95% of the dispersing ${\mathrm{SO}}_2$ plumes should be located within 300 km. Although this search is used to determine the wind-reconstruction, the stacked data spans over a box 1000 km × 1000 km. Over this 4-month eruptive period, SNPP and NOAA-20 together collected 493 plume snapshots comprising about $7.8\times {10}^6$ VCD estimates, generating the time-averaged plume in Figure 3 with fitting results summarized in

Table 1. Time-averaged quantities with uncertainty from fitting: 1 September–31 December 2018, Veniaminof volcano.

Parameter	Units	Median	[5th, 95th] Percentiles
${\mathrm{SO}}_2$ Flux	kt day${}^{-1}$	1.5	[0.8, 2.9]
wind speed	m s${}^{-1}$	12.0	[6.1, 23.4]
Eddy Diffusivity (downwind)	m${}^2$ s${}^{-1}$	8.2 × 10${}^5$	[4.2,16.0] × 10${}^5$
Eddy Diffusivity (crosswind)	m${}^2$ s${}^{-1}$	5.5 × 10${}^5$	[2.8, 10.8] × 10${}^5$
${\mathrm{SO}}_2$ lifetime	h	5.3	[2.7, 10.4]
Plume Peclet Number	-	3.4	[3.2, 3.5]

. Although the VCD values used here are estimated from the gridded CrIS ${\mathrm{SO}}_2$ z-scores (Figure 4a) and a conversion factor generated only from the most robust retrievals, there are enough samples to make a credible time-average. Despite this, there is still a subtle bend in the plume axis, which is the result of accumulating many plumes which themselves have been dispersed in complex wind fields (e.g., the individual plumes in Figure 4a). This crosswind asymmetry justifies the theoretical benefits of performing the fitting on an APD profile (Figure 3 inset) rather than a VCD grid.

3.2. Time Series of 10-Day Running Averages

Although this capability is useful for long-term ${\mathrm{SO}}_2$ flux estimation, using a shorter time-averaging window could make this estimation relevant to operational users such as volcano observatories. However, a shorter time-averaging window necessarily induces a lower signal-to-noise ratio. Because of the large volume of data used to make each average, even much shorter time-averages contain enough data to make usable estimates. Here we demonstrate the construction of a time series of 10-day averages for the same eruption from Veniaminof (Figure 4). The benefits of fitting an APD rather than VCD model are even more apparent in this context since each 10-day average VCD at a given location only contains about 40 CrIS VCD estimates, compared with the approximately 1500 CrIS estimates within 300 km of the plume axis.

During this eruptive period, there were principally four episodes of heightened ${\mathrm{SO}}_2$ emission. In some cases, these were characterized by a single short-duration, but sustained emission (late September), whereas other periods such as the increased ${\mathrm{SO}}_2$ flux in mid–late November were characterized by a sequence of pulses. In some cases, the peak fluxes do not coincide exactly with known periods explosions and ash emission, most notably, the significant explosive activity during 21–22 Nov. during which the AVO ACC/VAL were RED/Warning yet the ${\mathrm{SO}}_2$ flux peaked 3 days later. This example is most likely an artifact of the centered time-averaging (±5 days), which captured not only this significant ${\mathrm{SO}}_2$ release which itself lasted approximately 2 days, but also the fact that there were greater ${\mathrm{SO}}_2$ emissions following this event than preceding it. The largest ${\mathrm{SO}}_2$ emission occurred in mid-December with significant plumes dispersing north and northeast. This period of emission is discussed in greater detail below. Overall, the 4-month averaged flux (blue dashed line, Figure 4c) is substantially similar to the average of the time series of fluxes (blue dotted line, Figure 4c), although they are not exactly comparable owing to the statistics of lognormal variables. Although this capability yields intriguing details of ${\mathrm{SO}}_2$ emissions during this poorly understood eruption, a full analysis and geological interpretation of this eruptive sequence are beyond the scope of this work.

As described in Section 2.2.4, the ${\mathrm{SO}}_2$ decay rate and lifetime are approximately lognormal variables as predicted by theory and as evidenced by the distribution of ${\mathrm{SO}}_2$ lifetime estimates from a global catalogue of large ${\mathrm{SO}}_2$ sources [3]. At several points in the time series, particularly until mid-October, there are not enough substantial plumes in the time-averaging interval to make a credible estimate of the decay rate and lifetime (Figure 4e). For such intervals in the middle of the time series, the previous interval’s statistics (omitting flux) are used as defaults. For the earliest interval, we have used the global distribution of ${\mathrm{SO}}_2$ lifetimes from [3] to estimate the default lognormal parameters. As is evident later in the time series, the range of estimates for ${\mathrm{SO}}_2$ lifetime are very similar to the statistics of the lifetimes fitted by [3].

4. Discussion

4.1. Estimated Cloud Mass

Although geophysical parameters such as the flux, lifetime, and diffusivities are presented here as lognormal variables with uncertainty propagated directly from the decay rate estimation, the plume mass and Peclet number are direct result of the fitting and have uncertainty dominated by the fitting error. Owing to the statistics of lognormal variables, the median fitted plume mass can be calculated as the product of the median flux and median lifetime. Figure 4e shows the 10-day average cloud mass (blue) that would be estimated using only the detection steps introduced above, but not the stacking, that is, it is the 10-day average of the individual plume detections. This is a very conservative estimate since it does not include any ${\mathrm{SO}}_2$ below the z-score-based robust detection limit. For the plumes in this time series, this is about 40% of the fitted cloud mass (Figure 5).

Crucially, if the individually detected cloud masses could accurately incorporate regions of low-concentration ${\mathrm{SO}}_2$ surrounding the robust detections, then a plausible flux estimate could be made from a sequence of measurements simply by estimating the decay rate as in this work. However, clearly a significant amount of the total ${\mathrm{SO}}_2$ in these small clouds is too dilute for the robust detection and the rotation and stacking is required to amplify the signal in these regions.

4.2. Wind Speed, Goodness of Fit, and Anomalous Flux

For the 4 month-averaged plume, the fitted wind speed agrees very well with radiosonde wind data from Cold Bay, Alaska, located only 235 km to the southwest (Figure 6, data provided by Larry Oolman at the University of Wyoming Weather Lab, http://weather.uwyo.edu/upperair/sounding.html (accessed on 30 December 2020). Here we have analyzed all available soundings for the 4 month period, using only the data between 2 and 5 km altitude, which is the typical range of ${\mathrm{SO}}_2$ altitudes measured by CrIS during this period. The 4 month-averaged radiosonde wind speed was 13.0 m s${}^{-1}$ and the median (and 90% confidence interval) radiosonde wind was 11.8 (3.6, 26.8) m s${}^{-1}$ which compares very well with that from out fitted ${\mathrm{SO}}_2$ plume of 12.0 (6.1, 23.4) m s${}^{-1}$. At this large temporal scale, agreement is very good; however, this technique is not a robust tracer of time-averaged winds at finer scales.

As demonstrated by comparing Figure 4 and Figure 6b, the fitted wind speeds over 10-day intervals are generally significantly less than those from radiosonde winds. Principally, this discrepancy is the result of the fitting process based on ${\mathrm{SO}}_2$ plume shape, requiring a substantial enough plume to be present. Furthermore, the fitting is not a direct measure of the wind speed since any upwind ${\mathrm{SO}}_2$ dispersal (real or apparent and caused by smoothing) will greatly reduce the fitted wind speed. If more plume observations are incorporated, the time-averaged dispersal behavior emerges, yielding a much more accurate fitted wind speed. Fortunately, this discrepancy has little impact on the fitted ${\mathrm{SO}}_2$ flux and lifetime since lifetime is estimated independently and the flux is proportional to the fitted cloud mass, itself very accurate when the fit is good. When the fit is degraded, the fitted cloud mass (and thus flux) can be affected by the fitted wind speed.

Over the sequence of plumes from Veniaminof described here, the APD model fit is generally good with an $R^2\approx 0.9$ or better (Figure 4c). Critically, there is no apparent correlation between $R^2$ and flux; however, there is an apparent negative correlation between the wind speed estimate and $R^2$ for the most wind-affected plumes (Figure 4d). This typically results when the APD profile does not decrease smoothly in the downwind direction despite a very good fit (qualitatively) in the upwind direction. Owing to the simplified physics used, it could be possible to make a more robust estimate of the Peclet number and flux using only the physics of upwind dispersal; however, such a specialized fit is beyond the scope of this work.

As described, the largest ${\mathrm{SO}}_2$ emission occurred in mid-December, reaching peak emission rates of 8.6 kt day${}^{-1}$. Through this period, the fit worsened somewhat to $0.8<R^2<0.9$. Although many significant plumes can be seen in individual CrIS overpasses during this time (Figure 4a), it is possible that the fluxes during this period are overestimated as a result of an overestimated wind speed. If the wind has carried ${\mathrm{SO}}_2$ at an irregular speed or if the interval is characterized by a few large, discrete ${\mathrm{SO}}_2$ releases, it is very common to observe highly irregular (non-smooth) decay in the downwind direction. As the APD model attempts to fit the far-downwind ${\mathrm{SO}}_2$ that has not dispersed or diluted, the resulting fit is an apparently more wind-affected plume containing greater total mass and thus greater apparent flux with a poorer fit (Figure 7). Regardless of these considerations, the fitted model is still plausible throughout this period suggesting that the true fluxes were likely similar to those fitted here.

4.3. Comparison of Single-Sensor Results

Throughout, we have combined retrievals from NOAA-20 CrIS and SNPP CrIS. Here we demonstrate that these two sensors make consistent measurements by generating a 4-month averaged plume and fit for each sensor separately (Figure 8). The time-averaged VCD maps are largely similar with the NOAA-20 CrIS version extending over a slightly larger area downstream and the SNPP CrIS plume narrowing downstream and containing smaller peak VCD values than the NOAA-20 version. The fitted parameter values (

Table 2. Fitting results for NOAA-20 and SNPP-only time-averaging: 1 September–31 December 2018, Veniaminof volcano.

Parameter	Units	Percentiles: Median [5th, 95th]
		NOAA-20 Only	SNPP Only
${\mathrm{SO}}_2$ Flux	kt day${}^{-1}$	1.6 [0.8, 3.2]	1.0 [0.5, 2.0]
wind speed	m s${}^{-1}$	12.5 [6.4, 24.4]	8.1 [4.1, 15.8]
Eddy Diffusivity (downwind)	m${}^2$ s${}^{-1}$	7.2 [3.7, 14.0] × 10${}^5$	7.7 [4.0, 15.2] ×10${}^5$
Eddy Diffusivity (crosswind)	m${}^2$ s${}^{-1}$	2.6 [1.4, 5.2] × 10${}^5$	2.2 [1.1, 4.3] ×10${}^5$
${\mathrm{SO}}_2$ lifetime	h	5.1 [2.6, 10.0]	6.6 [3.4, 12.9]
Plume Peclet Number	-	4.0 [3.8, 4.2]	2.0 [1.9, 2.1]

) reflect these observations. The eddy diffusivities estimate for the combined analysis was significantly larger than for the individual sensors, in particular in the crosswind direction. This is due to the fact that combining the two sensors worth of observations can only increase the spread of the time-averaged VCD, mimicking the effect of mixing by large eddies. All other parameters were constrained between the NOAA-20 and SNPP fits; however, the combined fitted parameters are more similar to those of NOAA-20 than SNPP despite the fact that there were roughly the same number of NOAA-20 and SNPP observations.

4.4. Comparison with Other Estimates of Degassing and Lifetime

Although there are not any direct comparisons possible for this case study due to the limited coverage by UV instruments which typically make these measurments, we note here two promising comparisons that are available. In analyzing ${\mathrm{SO}}_2$ flux estimates from TROPOMI data, [13] compared annual emissions from a variety of large ${\mathrm{SO}}_2$ sources in the interval April 2018 to March 2019, including Veniaminof. They measured the annually averaged flux from Veniaminof in this period as $189\pm 14.7$ kt yr${}^{-1}$ (TROPOMI), $107\pm 27.5$ kt yr${}^{-1}$ (OMI), and $95\pm 31.6$ kt yr${}^{-1}$ (OMPS) which are in decreasing order of resolution and average sensitivity per area. Based on the 4-month CrIS-based estimate of $1.5$ kt day${}^{-1}$, CrIS measured a total of 183 kt ${\mathrm{SO}}_2$ over the course of the 2018 eruption. If very little or no emissions occurred in the other 8 months of the year, this would compare well with the estimate from TROPOMI in [13]. However, over the previous decade (2008–2017), annual average emissions from Veniaminof were approximately $62.1\pm 6.3$ kt yr${}^{-1}$ as measured by OMI for NASA’s public archive https://so2.gsfc.nasa.gov/measures.html (accessed on 29 December 2020) [38]. If this low level of emissions had occurred through the rest of the year, the annualized emissions from Veniaminof would be approximately $224\pm 83$ kt yr${}^{-1}$. Depending on which of the above UV measurements is most similar in sensitivity to CrIS, this annualized flux is between 19% and 236% greater that what was recorded with UV instruments although this comparison is probably most appropriate for CrIS-OMI (209%) and CrIS-OMPS (236%). Because about 100 kt of the ${\mathrm{SO}}_2$ was emitted in December when SZA was too high to make reliable UV measurements for all but the strongest plumes, it should be expected that CrIS would measure about double the annual ${\mathrm{SO}}_2$ flux as OMI and OMPS including a period of unrest occurring in the winter.

Lastly, as described throughout, the ${\mathrm{SO}}_2$ lifetime or decay time scale plays a central role in determining the ${\mathrm{SO}}_2$ flux. Consequently, our estimation technique for the decay rate outlined in Appendix E must be compared to other estimates. The most comprehensive data on ${\mathrm{SO}}_2$ plume parameters is compiled in the global catalog of large ${\mathrm{SO}}_2$ sources which includes a mix of volcanoes, coal-fired power plants, and smelting facilities [3]. The values of ${\mathrm{SO}}_2$ decay time scale in this catalog are the result of a similar wind-rotating and time-averaged fitting process, except that in [3], the wind speed is estimated independently of the fit using ancillary wind datasets and the decay time is fitted from the mean plume. Here, we adopt a completely different approach, estimating the ${\mathrm{SO}}_2$ lifetime from the record of individual plume mass observations, producing a lognormal probability distribution for the lifetime representing the uncertainty propagated from the noisy decay rate estimation process. Our distribution of ${\mathrm{SO}}_2$ lifetimes is in very good agreement with the distribution from the global catalog which contains data from 215 large ${\mathrm{SO}}_2$ sources (Figure 9). In particular, our distribution strays no more that 1% outside the 90% confidence interval for the distribution of the global catalog. We have generated this estimate using 10,000 empirical distribution functions (EDF) of simulated lifetime sample sets (each with 215 members) drawn from the global catalog histogram.

The strength of this decay time scale technique is that it can be estimated solely from a sequence of plume observations; however, the uncertainty in the technique is significant and propagates to the other fitted quantities unless additional data is used to constrain the problem. Overall, the good agreement between the distributions of ${\mathrm{SO}}_2$ lifetime estimates from Veniaminof and those from large sources globally is an indication that this technique can produce similar lifetimes as are measured by other means and thus that our flux estimates are at least as accurate as in similar techniques.

5. Conclusions

(i)

Time-averaged volcanic ${\mathrm{SO}}_2$ flux and lifetime estimation techniques can be successfully modified to accept newly-available hyperspectral IR ${\mathrm{SO}}_2$ data, enabling robust gas monitoring at high-latitude volcanoes in all illumination conditions, including low- or no sunlight conditions. This has the potential to fill a key observational gap in traditionally UV-based monitoring.

(ii)

Our technique detects weak signals of plumes in gridded ${\mathrm{SO}}_2$ data using computer vision and object-detection techniques used widely in medical and seismic tomography. It then amplifies these weak plume signals by attempting to reconstruct the wind fields dispersing them and time-averaging the results. This amplification enables the recovery of approximately twice as much ${\mathrm{SO}}_2$ mass as would be estimated by considering the cloud mass in individual plume snapshots alone. A new line-density fitting function derived directly from the time-averaged point-source dispersion physics is used to generate a well-fitting model of the time-averaged plumes.

(iii)

The ${\mathrm{SO}}_2$ decay rate and lifetime (reciprocal decay rate) are estimated independently of the fitting and are shown to have lognormal statistics which agree remarkably well with previous published estimates of the distribution of ${\mathrm{SO}}_2$ lifetimes from large sources globally.

(iv)

This technique can be subdivided, allowing for short-timescale averages and the generation of ${\mathrm{SO}}_2$ flux monitoring time series, which could prove useful to observatories. The technique in general and the time series capability in particular have been demonstrated for the September–December 2018 period of unrest and eruption at Veniaminof volcano, Alaska, highlighting the importance of new algorithms designed for detection and characterization of ${\mathrm{SO}}_2$ in the infrared.