Probabilistic retrieval of volcanic SO₂ layer height and partial column density using the Cross-track Infrared Sounder (CrIS)

During most volcanic eruptions and many periods of volcanic unrest, detectable quantities of sulfur dioxide (SO₂) are injected into the atmosphere at a wide range of altitudes, from ground level to the lower stratosphere. Because the fine ash fraction of a volcanic plume is, at times, colocated with SO₂ emissions, global tracking of volcanic SO₂ is useful in tracking the hazard long after ash detection becomes dominated by noise. Typically, retrievals of SO₂ vertical column density (VCD) have relied heavily on hyperspectral ultraviolet measurements. More recently, infrared sounders have provided additional VCD measurements and estimates of the SO₂ layer altitude, adding significant value to real-time monitoring of volcanic emissions and climatological analyses. These methods can provide fast and accurate physics-based retrievals of VCD and altitude without regard to solar irradiance, meaning that they are effective day and night and can observe high-latitude SO₂ even in the winter.

In this study, we detail a probabilistic enhancement of an infrared SO₂ retrieval method, based on a modified trace gas retrieval, to estimate SO₂ VCD and altitude probabilistically using the Cross-track Infrared Sounder (CrIS) on the Joint Polar Satellite System (JPSS) series of satellites. The methodology requires the characterization of real SO₂-free spectra aggregated seasonally and spatially. The probabilistic approach replaces altitude and VCD estimates with probability density functions for the layer height and the partial VCD at multiple heights, fully quantifying the retrieval uncertainty and allowing the estimation of SO₂ partitioning by layer. This framework adds significant value over basic VCD and altitude retrieval because it can be used to assign probabilities of SO₂ occurrence to different atmospheric intervals.

We highlight analyses of several recent significant eruptions, including the 22 June 2019 eruption of Raikoke volcano, in the Kuril Islands; the mid-December 2016 eruption of Bogoslof volcano, in the Aleutian Islands; and the 26 June 2018 eruption of Sierra Negra volcano, in the Galapagos Islands. This retrieval method is currently being implemented in the VOLcanic Cloud Analysis Toolkit (VOLCAT), where it will be used to generate additional cloud object properties for real-time detection, probabilistic characterization, and tracking of volcanic clouds in support of aviation safety.

1 Introduction

During most volcanic eruptions and many periods of volcanic unrest, detectable quantities of sulfur dioxide (SO₂) are injected into the atmosphere at a wide range of altitudes, from ground level to the lower stratosphere. Often early in eruptions the fine ash fraction of a volcanic plume is colocated with SO₂ emissions, and thus ash tracking can be performed by proxy; however, later ash and SO₂ tend to evolve along different trajectories due to subtly differences in altitude and removal processes (Karagulian et al., 2010; Corradini et al., 2010; Sears et al., 2013; Moxnes et al., 2014). Early colocation of SO₂ and ash is highly significant for informing forward trajectory models (e.g., HYSPLIT) of volcanic clouds as is performed in response to Volcanic Ash Advisories (VAAs) reported by the global network of Volcanic Ash Advisory Centers (VAACs). Because fine ash and SO₂ eventually diverge along different trajectories, due in large part to wind shear, layer height estimates are critical for ash and SO₂ cloud estimates. Although volcanic ash presents a demonstrated threat to aviation (e.g., ICAO, 2012; Casadevall, 1994; Prata and Rose, 2015; Guffanti et al., 2010), SO₂ also presents an aviation safety concern, mainly as a human health hazard and through damage by sulfuric acid, as well as via impacts on global climate and air quality (Chin and Jacob, 1996; Prata, 2009; Carn et al., 2009; Robock, 2000).

Globally, measurements of SO₂ vertical column density (VCD, here in Dobson Units, DU, 1 DU $=\mathrm{2.69}\times {\mathrm{10}}^{\mathrm{16}}$ molec. cm⁻²) have previously relied heavily on low-earth-orbiting hyperspectral ultraviolet (UV) instruments including the Ozone Monitoring Instrument (OMI), Ozone Mapping and Profiler Suite (OMPS), and Global Ozone Monitoring Experiment-2 (GOME-2) (e.g., Krotkov et al., 2010; Carn et al., 2017; Li et al., 2017; Theys et al., 2013). More recently, efforts to improve UV methods have focused on high-cadence UV measurements made from the Deep Space Climate Observatory – Earth Polychromatic Imaging Camera (DSCOVR-EPIC) at the Earth–Sun (L1) Lagrange point (Carn et al., 2018), as well as from high spatial resolution SO₂ VCD and limited layer height measurements from the Tropospheric Monitoring Instrument (TROPOMI) (Theys et al., 2019; Hedelt et al., 2019). In the last decade, infrared sounders such as the Infrared Atmospheric Sounding Interferometer (IASI) have provided additional SO₂ VCD measurements and estimates of the layer altitude, providing significant added value to real-time monitoring of volcanic emissions and climatological analyses (Walker et al., 2011, 2012; Carboni et al., 2012, 2016; Clarisse et al., 2014; Bauduin et al., 2016). These methods can provide fast and accurate physics-based retrievals of VCD and altitude. Furthermore, because these techniques principally rely on thermal contrast in the atmosphere and not solar irradiance (as in UV measurements), they are effective day and night and can observe high-latitude SO₂ in the winter months when UV techniques are unavailable. This twice-daily global coverage makes IR-based SO₂ retrieval a highly useful tool for operational support of aviation safety and for creating a truly continuous global analysis of SO₂ from volcanic eruptions.

In this study, we detail a probabilistic enhancement of the infrared SO₂ retrieval method of Clarisse et al. (2014), based on a modified trace gas retrieval (Walker et al., 2011) to estimate SO₂ VCD and altitude probabilistically utilizing the Cross-track Infrared Sounder (CrIS) currently aboard the Suomi-NPP (SNPP) and NOAA-20 satellites as part of the Joint Polar Satellite System (JPSS), having a local time ascending node (LTAN) of 13:30 with NOAA-20 operating approximately 50 min ahead of SNPP. Similar to IASI, CrIS is a Fourier transform Michelson interferometer covering three regions of the infrared spectrum: long-wave infrared (LWIR) (650–1095 cm⁻¹), mid-wave IR (MWIR) (1210–1750 cm⁻¹), and short-wave IR (SWIR) (2155–2550 cm⁻¹) (Han et al., 2013). Of the two principal SO₂ absorption features in the infrared (ν₁: 1000–1200 cm⁻¹; ν₃: 1300–1410 cm⁻¹), only the ν₃ band is covered by CrIS (MWIR). As highlighted by Carboni et al. (2012) and Clarisse et al. (2014), ν₃ is the stronger of the two absorption bands; however, it does contain significant interference from water vapor, limiting this retrieval's ability to characterize SO₂ features at very low altitudes. However, it is exactly the variable amounts of interference with water vapor at different heights that gives this technique its ability to retrieve SO₂ altitude information (Clarisse et al., 2014). As these studies pointed out, although clouds pose similar absorption features to water vapor, it is only in a broadband sense, allowing the finer absorption lines of SO₂ to be distinguished even in scenes with overlying meteorological clouds, as long as the clouds are not nearly opaque. Although the interference from water vapor is a significant theoretical limitation of this approach, in practice we have been able to extract some information on low-altitude SO₂ clouds even in the tropics which is detailed later. Despite these limitations, the ν₃ band absorption lines are only minimally influenced by ash and dust, making this height retrieval method especially useful early in the evolution of volcanic eruption clouds where there is typically colocation between SO₂ and ash clouds.

Both CrIS instruments are currently operating in full spectral resolution mode (FSR), providing MWIR spectra at 0.625 cm⁻¹ spectral resolution since December 2015 for SNPP CrIS (excepting a major outage 26 March–1 August 2019) and February 2018 for NOAA-20 CrIS. CrIS scans consist of 30 fields of regard (FOR) in 3.3^∘ steps between ±48.3^∘ scan angle, each of which contains 9 circular fields of view (FOV) arranged in a square (3×3) array that rotates and stretches as the mirror moves away from nadir towards edge of scan (Han et al., 2013). CrIS granules are collected into 6 min granules of 45 scans, resulting in 12 150 MWIR FSR spectra collected every 6 min. The CrIS swath width is 2200 km; however, because of the rotating FORs, some ground points are measured by multiple FOVs even within the same scan, and some gaps exist due to the square FOR layout of circular FOVs and the presence of short gaps between scans. The FOV at the center of each FOR (number 5) is a 14 km diameter circle at nadir, extending out to an 43.6 km × 23.2 km (major and minor axes) ellipse for the first and last FORs on the edge of the swath (Han et al., 2013; Wang et al., 2013). Although CrIS FOVs are slightly larger than IASI FOVs (14 km versus 12 km at nadir), there are many more of them per scan since CrIS FOR are 3×3 arrays, whereas IASI FOR are 2×2 with larger gaps between FOV, FORs, and scan lines but the same swath width (e.g., Sun et al., 2018). Consequently, CrIS makes many more measurements per area than IASI does, resulting in greater overall resolution than IASI. Lastly, all of the CrIS MWIR channels used in the present study are, in general, very low noise (noise equivalent differential radiance, NEdN <0.05 mW m⁻² sr⁻¹ cm) with the exception of SNPP CrIS FOV 7, which is above specification. Later in this study we will show some retrievals from this FOV; however, these are considered to be of very low quality and are not considered reliable. They are shown here only to elucidate how strong instrument noise is propagated to the retrieval.

Figure 1 Flowchart showing the probabilistic framework for Monte Carlo height and VCD, yielding a PDF for the height, which is not generally Gaussian and may be heavily skewed, and a Gaussian distribution of conditional VCD $\left(\widehat{X}\right|H=h)$. The height retrieval of Clarisse et al. (2014) is shown schematically as red lines, giving a single height estimate, which is not the mean height in general (approximately the dotted black line in the height PDF). As this figure is only a schematic, the pictorial relationship between these height estimates in not universal.

The NOAA Unique CrIS/ATMS Processing System (NUCAPS) already includes a retrieval of SO₂ from CrIS data (Gambacorta, 2013); however, it is based on a heritage algorithm designed to estimate many trace gases from cloud-cleared radiances in one retrieval, whereas we focus more specifically on the problem of retrieving SO₂ in any background atmosphere from all available CrIS measurements. The methodology requires the characterization of the background mid-wave infrared spectrum of the SO₂-free atmosphere, which is done by collecting the statistics of more than 360 million SO₂-free CrIS spectra aggregated seasonally and spatially. The probabilistic approach replaces altitude and VCD estimates with a nonparametric probability density function (PDF) for the layer height and estimates (with uncertainty) of the partial VCD at multiple heights, fully quantifying the retrieval uncertainty and allowing the estimation of SO₂ partitioning by layer (Fig. 1). This framework adds significant value because it can be used to assign probabilities of SO₂ occurrence in different intervals of the atmosphere, which could prove very useful for aviation safety in the future when changing aviation hazard priorities will require such information (ICAO, 2019).

In this study we analyze several recent significant eruptions, including the 22 June 2019 eruption of Raikoke volcano, in the Kuril Islands; the mid-December 2016 eruption of Bogoslof volcano, in the Aleutian Islands; and the 26 June 2018 eruption of Sierra Negra volcano, in the Galapagos Islands. This retrieval method is currently being implemented in the VOLcanic Cloud Analysis Toolkit (VOLCAT, https://volcano.ssec.wisc.edu/, last access: 26 October 2020; Pavolonis et al., 2013, 2015a, b, 2018), where it will be used to generate additional cloud object properties for real-time detection, probabilistic characterization, and tracking of volcanic clouds in support of aviation safety.

2 Probabilistic SO₂ layer retrieval theory

2.1 Classical methods for height retrieval

As a preliminary we discuss several methods which we describe here as “classical”. In fact these methods are relatively recent; however, they do not make full use of the probability spaces that we will exploit here. Previous analyses of the height and distribution of volcanic SO₂ plumes using data from IASI by Carboni et al. (2012, 2016) and Clarisse et al. (2014) utilized trace gas methods modified from the method originally outlined by Walker et al. (2011). The analysis of Carboni et al. (2012) imposed an a priori Gaussian vertical distribution over pressure coordinates for the SO₂ concentration, retrieving the total SO₂ VCD, the mean pressure, and the standard deviation pressure (spread, only if VCD sufficiently strong). By contrast, Clarisse et al. (2014) developed a system in which the SO₂ is assumed to exist in a narrow box profile layer and an iterative retrieval is performed for the VCD concentration conditional on the retrieved SO₂ layer altitude. The principal differences between these methods and the method detailed here are summarized in Table 1. Employing the notation of Rodgers (2000), this retrieval relates a set of parameters governing the concentration of a trace gas (the true state, x∈ℝ^M, SO₂ in this case) to a set of measurements y∈ℝ^N (typically brightness temperature spectra) by a forward radiative transfer model $\mathit{F}:{\mathbb{R}}^M\to {\mathbb{R}}^N$. In the following exposition, we focus on and enhance the method of Clarisse et al. (2014).

Carboni et al. (2012)Clarisse et al. (2014)

Table 1 Summary of Recent Infrared SO₂ Height Estimation Methods.

The infrared trace gas methods of Walker et al. (2011) and Clarisse et al. (2014) rely on the ability to write the data and true state each as a sum of their respective climatological background averages (${\stackrel{\mathrm{‾}}{\mathit{y}}}_{\mathrm{bg}},{\stackrel{\mathrm{‾}}{\mathit{x}}}_{\mathrm{bg}}$) and their anomalies ($\stackrel{\mathrm{̃}}{\mathit{y}},\stackrel{\mathrm{̃}}{\mathit{x}}$). Linearizing around the climatological average gives

$$\begin{array}{}\text{(1)}& {\stackrel{\mathrm{‾}}{\mathit{y}}}_{\mathrm{bg}}+\stackrel{\mathrm{̃}}{\mathit{y}}=\mathit{F}\left({\stackrel{\mathrm{‾}}{\mathit{x}}}_{\mathrm{bg}};\mathit{u}\right)+\mathbf{K}\stackrel{\mathrm{̃}}{\mathit{x}}+{\mathit{\epsilon }}_{\mathrm{tot}},\end{array}$$

where u is a collection of all auxiliary parameters necessary to the forward model including atmospheric pressure, temperature, and water vapor profiles, as well as the state of the surface and instrument. $\mathbf{K}\in {\mathbb{R}}^{N\times M}$ is the Jacobian of F linearized around ${\stackrel{\mathrm{‾}}{\mathit{x}}}_{\mathrm{bg}}$, and ε_tot is the total error associated with the measurement, linearization, and other surface and atmospheric properties (including clouds) that influence the measured radiances. Because it is inferred that ${\stackrel{\mathrm{‾}}{\mathit{y}}}_{\mathrm{bg}}=\mathit{F}({\stackrel{\mathrm{‾}}{\mathit{x}}}_{\mathrm{bg}};\mathit{u})$, the equation for the error is reduced simply to a relationship between the state and measurement anomalies:

$$\begin{array}{}\text{(2)}& {\mathit{\epsilon }}_{\mathrm{tot}}=\stackrel{\mathrm{̃}}{\mathit{y}}-\mathbf{K}\stackrel{\mathrm{̃}}{\mathit{x}}.\end{array}$$

From this formula, the retrieval of $\stackrel{\mathrm{̃}}{\mathit{x}}$ can proceed either by maximum likelihood estimation, iterative methods such as Levenberg–Marquardt and gradient descent algorithms, or other methods.

For a 1 km thick box profile layer, the concentration of anomalous SO₂ can be represented by two parameters: the total VCD of the gas (x) and the height of the layer center (h). Using such a profile, the model spectrum is then a function of the VCD and the layer height. In Clarisse et al. (2014) and this paper, the Jacobian is represented a function of height but only contains a VCD perturbation. Because this Jacobian only measures the sensitivity of the forward model to the presence of SO₂ at each height independently, it is best viewed as a set of vectors in the same space as y rather than as a matrix and thus we write it hereafter as K(h)∈ℝ^N. Additionally, the model Jacobian for the trace gas retrieval is calculated at the background state (zero for SO₂), and the Jacobian is approximated as a finite difference at each altitude:

$$\begin{array}{}\text{(3)}& \mathit{K}\left(h\right)\approx {\displaystyle \frac{\mathit{F}(\mathit{ϵ},h;\mathit{u})-\mathit{F}(\mathrm{0},h;\mathit{u})}{\mathit{ϵ}}},\end{array}$$

where the perturbation (ϵ) is taken as 5 DU.

Figure 2 (a) Channel-wise brightness temperature difference response to increasing VCD. All channels here are used in the original retrieval. Red channels are used in the specialized retrieval due to their mostly linear response across the full range of reasonable VCD. (b) Height-dependent Jacobians (normalized) showing specialized channel selection in panel (a).

In the present study, we pre-compute a limited database of Jacobians for 1 km thick SO₂ layers centered at 28 altitudes between 0 and 32 km (Fig. 2b; Fig. 1 of Clarisse et al., 2014). We use standard profiles of pressure, temperature, and water vapor for a tropical atmosphere, summer and winter midlatitude atmospheres, and summer and winter subarctic atmospheres, resulting in a total of 140 Jacobians to be used in the retrieval (Clough et al., 2005). As in Pavolonis (2010), all radiative transfer model simulations used here were performed using the LBLDIS tool (Turner, 2005), which utilizes the Line-by-Line Radiative Transfer Model (LBLRTM; Clough and Iacono, 1995; Clough et al., 2005) to compute gaseous absorption and the Discrete Ordinate Radiative Transfer (DISORT) model to complete the radiative transfer calculation (including multiple scattering).

Following Clarisse et al. (2014), a height-dependent Jacobian can be used to calculate a statistical z score, measuring the relative confidence in the presence of SO₂ at each height:

$$\begin{array}{}\text{(4)}& z\left(h;\mathit{y},{\stackrel{\mathrm{‾}}{\mathit{y}}}_{\mathrm{bg}}\right)={\left[{\mathit{K}}^{\mathrm{T}}\left(h\right){\mathbf{S}}^{-\mathrm{1}}\mathit{K}\left(h\right)\right]}^{-\frac{\mathrm{1}}{\mathrm{2}}}{\mathit{K}}^{\mathrm{T}}\left(h\right){\mathbf{S}}^{-\mathrm{1}}\left(\mathit{y}-{\stackrel{\mathrm{‾}}{\mathit{y}}}_{\mathrm{bg}}\right),\end{array}$$

where the mean climatological background spectrum ${\stackrel{\mathrm{‾}}{\mathit{y}}}_{\mathrm{bg}}$ and error covariance matrix S are built from spectral residuals from a database of measurements with low or no detectable SO₂ present. The construction of this database of SO₂-free spectra is detailed in Sect. 2.5. Note that because K(h) is a vector, the factor K^T(h)S⁻¹K(h) is a scalar at each height h. Here, $z(h;\mathit{y},{\stackrel{\mathrm{‾}}{\mathit{y}}}_{\mathrm{bg}})$ is the statistical z score (number of standard deviations from the mean) of finding the SO₂ anomaly at altitude h given the data y and the SO₂-free background spectrum ${\stackrel{\mathrm{‾}}{\mathit{y}}}_{\mathrm{bg}}$. Using the z score, Clarisse et al. (2014) estimated the layer height (which we refer to as h_C) as that which maximizes the z-score function:

$$\begin{array}{}\text{(5)}& h_{\mathrm{C}}:=\underset{h}{\arg \max }z\left(h;\mathit{y},{\stackrel{\mathrm{‾}}{\mathit{y}}}_{\mathrm{bg}}\right).\end{array}$$

This method was used in that study to produce a consistent and reasonably accurate set of cloud top height estimates for the 2011 eruption of Nabro Volcano, Eritrea. This method is currently the principal operational SO₂ layer height method for IASI used by the Support to Aviation Control Service (SACS, https://sacs.aeronomie.be/, last access: 26 October 2020) of the Royal Belgian Institute for Space Aeronomy (BIRA-IASB), supporting several Volcanic Ash Advisory Centers (VAACs) in near-real time (Brenot et al., 2014).

For simplicity, throughout the remainder of this work we refer to this type of height retrieval with a function notation:

$$\begin{array}{}\text{(6)}& h_{\mathrm{C}}=g\left(\mathit{y},{\stackrel{\mathrm{‾}}{\mathit{y}}}_{\mathrm{bg}}\right):=\underset{h}{\arg \max }z\left(h;\mathit{y},{\stackrel{\mathrm{‾}}{\mathit{y}}}_{\mathrm{bg}}\right).\end{array}$$

2.2 Probabilistic enhancement

Throughout, we make a distinction between our method being probabilistic and other methods being deterministic; however, we note here that the classical methods are all based on the optimal estimation (a Bayesian method) of Rodgers (2000) and therefore are probabilistic in the sense that the retrieved quantities are the (Gaussian) mean and covariance of the state estimate. We make this distinction to highlight the fact that in the present study we focus on a much more detailed uncertainty propagation, specifically, propagation uncertainty about the SO₂-free background atmosphere to the height retrieval, allowing some departure from the Gaussian assumption underlying previous methods. Although the Gaussian assumption is workable for many types of retrieval, it is unsuitable for the Clarisse et al. (2014) height retrieval in particular due to the role played by the arg max operation. To see this directly, we must “probabilize” the Clarisse et al. (2014) height retrieval as follows. In this process, we use a notation common in probability theory in which random variables are represented as capitalized versions of their deterministic realizations. In what follows, the only exception to this notation will be that the forward model, Jacobian, and covariance matrix are not random variables.

Instead of calculating the mean and covariance of the climatological background as in a traditional trace gas retrieval, we treat the background (SO₂-free) spectrum as a random vector Y_bg (rather than the realization ${\stackrel{\mathrm{‾}}{\mathit{y}}}_{\mathrm{bg}}$, which is its mean), where the vector elements (each of the sampled wavenumbers or channels) are random variables. Although the total uncertainty in the measured SO₂-free background spectrum contains a mix of aleatoric (stochastic in fact) and epistemic (knowledge-deficit) uncertainties, the epistemic uncertainty due to a lack of knowledge about the true SO₂-free background for a given spectrum is considerably larger than aleatoric sources such as instrument noise and scattering. Consequently, we generate a probability space for the SO₂-free spectrum Y_bg, the uncertainty for which is derived principally from our ignorance of the true atmospheric state if SO₂ were not present. The brightness temperatures for each channel (element of Y_bg) is characterized by its own marginal probability distribution. In reality, these distributions may belong to a family of parameterized distributions; however, in this study they are nonparametric and only characterized by a marginal distribution (histogram) on each channel. Additionally, the elements of Y_bg are correlated, which is given by the covariance structure:

$$\begin{array}{}\text{(7)}& \mathbf{S}=\mathbb{E}\left[\left({\mathit{Y}}_{\mathrm{bg}}-\mathbb{E}\left[{\mathit{Y}}_{\mathrm{bg}}\right]\right){\left(\mathit{Y}-\mathbb{E}\left[{\mathit{Y}}_{\mathrm{bg}}\right]\right)}^{\mathrm{T}}\right],\end{array}$$

where the expectation is an average over all SO₂-free spectra in the database (detailed in Sect. 2.5).

In this framework, the z-score function is a conditional random variable given the layer height h:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}Z& \left(h;\mathit{y},{\mathit{Y}}_{\mathrm{bg}}\right)=Z|H=h\\ & ={\left[{\mathit{K}}^{\mathrm{T}}\left(h\right){\mathbf{S}}^{-\mathrm{1}}\mathit{K}\left(h\right)\right]}^{-\frac{\mathrm{1}}{\mathrm{2}}}{\mathit{K}}^{\mathrm{T}}\left(h\right){\mathbf{S}}^{-\mathrm{1}}\left(\mathit{y}-{\mathit{Y}}_{\mathrm{bg}}\right),\end{array}\end{array}$$

and the height is therefore a random variable $H=g(\mathit{y},{\mathit{Y}}_{\mathrm{bg}})$.

Implicitly, Clarisse et al. (2014) assumes that Y_bg is a multivariate normal random vector with mean ${\stackrel{\mathrm{‾}}{\mathit{y}}}_{\mathrm{bg}}=\mathbb{E}\left[{\mathit{Y}}_{\mathrm{bg}}\right]$ and covariance S, meaning that Z is a standard normal random variable. This fact about the z score is expected to hold in the present case with the full probabilistic characterization of the generally non-Gaussian Y_bg because the z score is a weighted sum over all of the channels in Y_bg, which is expected to converge to a Gaussian for a large collection of channels. Because the function g uses the arg max operation, which in not exactly a proper function (and also not linear), we can write that

$$\begin{array}{}\text{(9)}& \mathbb{E}\left[g\left(\mathit{y},{\mathit{Y}}_{\mathrm{bg}}\right)\right]\ne g\left(\mathit{y},\mathbb{E}\left[{\mathit{Y}}_{\mathrm{bg}}\right]\right).\end{array}$$

That is, the random variable resulting from a nonlinear transformation of a Gaussian random variable is not itself Gaussian and the mean value of that new variable is not equal to the value obtained by transforming the mean of the Gaussian (thus, 𝔼[H]≠h_C; Fig. 1 schematically), a standard result in elementary probability theory texts (e.g., DeGroot and Schervish, 2012). Similarly, h_C is not generally the maximum likelihood height either (h_C≠mode[H]). Consequently, without a clear understanding of what h_C is measuring in terms of the statistics of H, it is difficult to contextualize the value h_C. The principal enhancement over the classical method comes from setting the height retrieval in a probabilistic framework, enabling precise propagation of uncertainty in the background state to uncertainty in the retrieved height.

This study aims to estimate the probability distribution of H and show the importance of its PDF in making predictions about the cloud. We enhance the method of Clarisse et al. (2014) by retrieving a probabilistic SO₂ layer; i.e., we retrieve the SO₂ layer height as a PDF for the height (Fig. 1). As described, the probabilistic nature of the retrieval product is derived from propagating our uncertainty about the SO₂-free background spectrum through the Clarisse et al. (2014) height retrieval method (Fig. 1). Here, we use a set of 10 000 possible SO₂-free background spectra computed by Monte Carlo (MC) sampling according to the collection of marginal distributions of Y_bg and its covariance matrix S, which are first computed from a database of SO₂-free spectra (detailed in Sect. 2.5). The process for sampling this generally non-Gaussian correlated random vector is detailed in Appendix A. In this study each sample is denoted as ${\mathit{y}}_{\mathrm{bg}}^s\in {\mathrm{\Omega }}_{{\mathit{Y}}_{\mathrm{bg}}}$, where ${\mathrm{\Omega }}_{{\mathit{Y}}_{\mathrm{bg}}}$ is the sample space of Y_bg.

Although we could directly estimate the height PDF from sampling the many different backgrounds, we treat our retrieval as an update on the Clarisse et al. (2014) height estimate and cast this process in a Bayesian framework. In this framework, we treat the estimate and uncertainty from the (Clarisse et al., 2014) method as a prior distribution for the height and construct an approximate likelihood function from the data not accounted for directly in that method (the distribution of the many possible spectral residuals). The height PDF which is sought is the posterior distribution.

We impose a Gaussian prior with mean and variance given by the Clarisse et al. (2014) method. To estimate the mean and variance of the Gaussian, we first retrieve the Clarisse et al. (2014) height h_C and then generate many height estimates around h_C using a model spectral anomaly with SO₂ assumed at h_C and MC sampling of the zero-mean noise contained in the collection of possible backgrounds. Specifically, we estimate the height due to noisy model spectral anomaly samples

$$\begin{array}{}\text{(10)}& {\stackrel{\mathrm{̃}}{\mathit{y}}}^s=\left[\mathit{F}\left(\mathit{ϵ},h_{\mathrm{C}},\mathit{u}\right)-\mathit{F}\left(\mathrm{0},h_{\mathrm{C}},\mathit{u}\right)\right]-\left[{\mathit{y}}_{\mathrm{bg}}^s-{\stackrel{\mathrm{‾}}{\mathit{y}}}_{\mathrm{bg}}\right].\end{array}$$

This anomaly represents a modeled spectral anomaly with zero-mean spectral noise added and allows for the possibility of a bias induced by the difference between the mean climatological background and the model background which does not include cloud layers. The Gaussian prior mean and variance are then taken to be the mean and variance of these noisy modeled samples. Of note, this Gaussian prior mean is very close to the value h_C, but is preferred since its use does not restrict the Gaussian to be centered only at height values for which Jacobians were computed. Using these values for mean and variance is very similar to the estimate with uncertainty shown in Fig. 2 of Clarisse et al. (2014). This mean and variance parameterizes the Gaussian prior distribution $f_H^{\mathrm{prior}}\left(h\right)$.

The likelihood function is constructed directly by retrieving the height due to the real spectrum y and the set of sampled background spectra. Each height sample is generated as follows:

$$\begin{array}{}\text{(11)}& h^s=g\left(\mathit{y},{\mathit{y}}_{\mathrm{bg}}^s\right),\end{array}$$

i.e., we construct the random variable from elements of its sample space h^s∈Ω_H. The likelihood function measures the distribution of the possible real spectral residuals (y−Y_bg) given an SO₂ layer at height h. As there is no analytic probability model to describe this, we estimate that the likelihood is proportional to the distribution of possible heights as computed by kernel density estimation (KDE, e.g., Silverman, 1986) on this set of height samples:

$$\begin{array}{}\text{(12)}& \mathcal{L}\left(h;\mathit{y}-{\mathit{Y}}_{\mathrm{bg}}\right):=\widehat{f}\left(\mathit{y}-{\mathit{Y}}_{\mathrm{bg}}|h\right)\propto \mathrm{KDE}\left(\mathit{\left\{}{h}^s\mathit{\right\}}\right).\end{array}$$

An alternative approach would be to replace the KDE distribution with a simple histogram of the samples, although we use KDE because we seek a distribution which is at least piecewise-continuous and not piecewise-constant. This gives an estimate of the posterior height PDF:

$$\begin{array}{}\text{(13)}& f_H\left(h\right)\propto \mathcal{L}\left(h;\mathit{y}-{\mathit{Y}}_{\mathrm{bg}}\right)f_H^{\mathrm{prior}}\left(h\right),\end{array}$$

where the proportionality is eliminated by normalizing the posterior PDF such that the total probability is unity.

Although slower than retrieving the layer height (h_C) due to a mean spectrum alone, this distribution provides significantly more information, including the full PDF of H. This PDF may be used to calculate the modal, mean, and median values of the retrieved height or probabilities of finding the plume in a given altitude interval. Additionally, this PDF is essential for calculating the VCD correctly according to probability theory as detailed in the following section.

2.3 Probabilistic vertical column density

Although this method is used primarily for detection (using z scores) and height estimation, we estimate VCD as a side-product, which we treat here as a random variable $\widehat{X}$. Because we will use a linearized retrieval, the VCD values presented here will not be as accurate as those produced by an iterative technique. However, the details of our process below produce VCD values that are reasonably accurate for all but very strong emissions, as is demonstrated later in this work. Specifically, the way in which the height PDF is incorporated (detailed below) mitigates some underestimation error that would otherwise occur in a linearized approach for even dilute SO₂ clouds. As with the height estimation, the uncertainty propagated to the VCD primarily represents uncertainty about the SO₂-free background, which is of great importance in these linearized trace gas methods.

Because the estimated VCD depends strongly on the layer height, we refer to an estimate of total VCD where the layer is given as a specified height as a “conditional VCD”. In this framework, the VCD estimates of Walker et al. (2011, 2012) are conditional VCDs and are represented as a conditional random variable:

$$\begin{array}{}\text{(14)}& \begin{array}{rl}\widehat{X}|H& =h\\ & =\cos \mathit{\theta }{\left[{\mathit{K}}^{\mathrm{T}}\left(h\right){\mathbf{S}}^{-\mathrm{1}}\mathit{K}\left(h\right)\right]}^{-\mathrm{1}}{\mathit{K}}^{\mathrm{T}}\left(h\right){\mathbf{S}}^{-\mathrm{1}}\left(\mathit{y}-{\mathit{Y}}_{\mathrm{bg}}\right),\end{array}\end{array}$$

where an air mass factor equal to the cosine of the satellite zenith angle (cos θ) has been applied. This formula is in some sense a probabilistic enhancement of an optimal unconstrained least-squares estimate. Similar to the z scores, this function is normally distributed at every height with mean $\mathbb{E}\left[\widehat{X}\right|H=h]$ and variance $\mathrm{Var}\left[\widehat{X}\right|H=h]$. These are calculated as the sample mean and variance of the conditional VCD samples due to the many possible background spectra used to estimate the height PDF. The conditional VCD is a random function of height which generally reflects the principles of the Beer–Lambert law for any given realization; i.e., a smaller VCD is retrieved for a given spectral anomaly if the layer is assumed to be higher in the atmosphere. If the height of the SO₂ layer were known exactly, the VCD could be estimated by evaluating the conditional VCD function at that exact height. This is exactly what is done in the VCD retrieval of Clarisse et al. (2014), except by using an iterative conditional VCD calculation. However, in this study, since the height of the layer is known only probabilistically, i.e., as measured by the PDF f_H(h), additional computation is required to determine the VCD.

Although we do not know the true vertical profile of SO₂ concentration, the retrieval assumes a thin layer representation of the SO₂. The total VCD ($\widehat{X}$, a random variable) is obtained by integrating the box profile between the ground and the top of the atmosphere. Similarly, a partial VCD, denoted here as $\widehat{X}\left(h\right)$, can be calculated by integrating the box profile between the ground and the height h and is zero for h<H and rises linearly within the layer to $\widehat{X}$ for h≥H. Because the assumed concentration profile scales linearly with the total VCD ($\widehat{X}$) and the conditional VCD is normally distributed at each height, the partial VCD is normally distributed as well, thus requiring only two parameters: the mean and variance, which can be found using the conditional VCD function calculated above.

Here, we give approximation formulae for the mean and variance partial VCD. The derivation of these formulae is detailed in Appendix B. The mean partial VCD below height h is found by the law of total expectation:

$$\begin{array}{}\text{(15)}& {\mathit{\mu }}_{\widehat{X}}\left(h\right):=\mathbb{E}\left[\widehat{X}\left(h\right)\right]=\underset{\mathrm{0}}{\overset{h}{\int }}{f}_H\left(\mathit{\eta }\right)\mathbb{E}\left[\widehat{X}|H=\mathit{\eta }\right]\mathrm{d}\mathit{\eta },\end{array}$$

where the expectation of the conditional VCD is taken as the mean of the samples. This formula represents a weighted average of mean conditional VCD values where the height probability density assigns the weights. Because the conditional VCD is a function of height, this formula mixes the VCD estimates from different assumed heights.

The variance can also be calculated from the statistics of the conditional VCD expectation:

$$\begin{array}{}\text{(16)}& \begin{array}{rl}{\mathit{\sigma }}_{\widehat{X}}^{\mathrm{2}}\left(h\right)& :=\mathrm{Var}\left[\widehat{X}\left(h\right)\right]\\ & =\underset{\mathrm{0}}{\overset{h}{\int }}{f}_H\left(\mathit{\eta }\right)[\mathrm{Var}\left[\widehat{X}|H=\mathit{\eta }\right]\\ & +{\left(\mathbb{E}\left[\widehat{X}|H=\mathit{\eta }\right]\right)}^{\mathrm{2}}]\mathrm{d}\mathit{\eta }-{\mathit{\mu }}_{\widehat{X}}^{\mathrm{2}}(h),\end{array}\end{array}$$

where the conditional mean and conditional variance were previously estimated from the MC samples.

The covariance between the partial VCDs for two altitudes (a and b) is given in Appendix B. Of particular interest is that these formulae may be used to calculate the expectation and variance values of the partial VCD between two altitudes:

$$\begin{array}{}\text{(17a)}& {\displaystyle }\mathbb{E}\left[\widehat{X}\left(b\right)-\widehat{X}\left(a\right)\right]={\mathit{\mu }}_{\widehat{X}}\left(b\right)-{\mathit{\mu }}_{\widehat{X}}\left(a\right),\text{(17b)}& {\displaystyle }\begin{array}{rl}\mathrm{Var}\left[\widehat{X}\left(b\right)-\widehat{X}\left(a\right)\right]& ={\mathit{\sigma }}_{\widehat{X}}^{\mathrm{2}}\left(b\right)\\ & +{\mathit{\mu }}_{\widehat{X}}\left(a\right)\left({\mathit{\mu }}_{\widehat{X}}\left(b\right)-{\mathit{\mu }}_{\widehat{X}}\left(a\right)\right).\end{array}\end{array}$$

In this system, we retrieve probabilistic SO₂ information in two stages. In the first stage we perform an initial detection using the classical method (Eq. 4) to pre-screen each CrIS FOV that likely contains SO₂, taken as an initial maximum z score greater than 5, i.e, $z(h_{\mathrm{C}};\mathit{y},{\stackrel{\mathrm{‾}}{\mathit{y}}}_{\mathrm{bg}})>\mathrm{5}$ (e.g., Walker et al., 2011, 2012; Clarisse et al., 2014). Preliminary investigation of this threshold indicates that it somewhat conservative, striking a balance between including regions of diffuse SO₂ and excluding almost all false detections. In the second stage, we retrieve the height PDF and the mean and variance partial VCD as a function of height for each CrIS FOV that satisfies this initial z-score threshold.

2.4 Specialized retrieval for strong SO₂ loading

For strong SO₂ columns, an alternate retrieval is needed to increase sensitivity of the retrieved VCD owing to error induced by the linearized retrieval. We define strong loading here heuristically as z>200, which in preliminary testing corresponded to VCD values between approximately 10 and 20 DU. Since the conditional VCD retrieval uses a linearized forward model with only a 5 DU perturbation, it is expected that such an approximation would only hold for values near 5 DU and that many physically realistic VCD values would fall outside the linearization's radius of convergence. For large VCD values, the sensitivity of linear Jacobians to additional SO₂ is greatly reduced for most CrIS channels, especially the strongest CrIS channels, so linearized Jacobians with only a 5 DU anomaly will drastically underpredict the VCD for a given brightness temperature difference (Fig. 2a). In order to construct a linearized Jacobian that is more sensitive at higher VCD values, two approaches are possible: (i) we could use a larger VCD perturbation (a coarser finite difference), or (ii) we could seek a special combination of channels for which the linearization is a good approximation. The first approach is limited in that the strongest responses in the forward model are highly nonlinear, so a coarse finite difference will induce large errors there. The second approach is promising as long as such a suitable channel selection can be made that still preserves some of the main features of the ν₃ absorption band, retains enough channels to be robust, and is only applied when it is certain that SO₂ is dominating the signal.

In addressing this issue, we adopt the second approach. The specialized Jacobian must be dominated by channels with approximately linear forward model responses (Fig. 2b). This was accomplished in practice by constructing a sequence of Jacobians with various finite difference coarseness and choosing those channels for which the sequence of Jacobians is approximately constant. This channel subset (Appendix D) is used with the original 5 DU Jacobian, which can then be extrapolated successfully to high VCD values because the forward model truly is approximately linear for those channels.

One complication here is the fact that the channels with the most linear response are also those which are least sensitive to SO₂ VCD (Fig. 2). However, by applying this new retrieval only when strong SO₂ loading (z>200) is detected by the pre-screening retrieval, the signal is guaranteed to be dominated by SO₂ absorption even in the weak, approximately linear channels, and the linearization is expected to have moderately good accuracy. In theory, a sequence of increasingly restricted retrievals could be implemented to increase the sensitivity to even stronger SO₂ loads; however, even the most sensitive channels in the specialized subset show the worst linear approximations (Fig. 2a), and so some underestimation in the inversion will always be present with such an approach. As we will demonstrate below, this approach does not increase the sensitivity so much that extremely high VCD values (several hundred to thousands of Dobson units) can be retrieved but instead increases the ability of the linearized approach to resolve moderate to high VCD values (perhaps tens to low hundreds of Dobson units).

2.5 Background state construction

At every stage of the retrieval, the background state of the volcanic SO₂-free atmosphere must be accurate in order for this linearized method to succeed. Consequently, calculating accurate statistics of the background state is paramount. Each CrIS instrument collects almost 3 million spectra per day, allowing for robust characterization of the background spectrum including variation in conditions across seasons and locations. Because we use real CrIS spectra, not only are meteorological variations (including water vapor, clouds, temperature, etc.) accounted for, but the instrument noise profile is also included.

Figure 3 (a) IASI-derived days with (blue) and without (white) SO₂ columns with VCD≥1 DU in the 1-year background construction interval (1 November 2017–1 November 2018). Date intervals across the top define the seasons. (b) Histogram showing the number of SNPP CrIS spectra in each latitude–longitude cell totaled over the 1-year interval.

In constructing the background spectrum channel-wise marginal PDFs (histograms) and covariance matrix, periods with little or no SO₂ must be determined. We utilize the detailed record of global volcanic SO₂ emissions from the operational IASI SO₂ retrieval algorithm (Lieven Clarisse, personal communication, 2018) between 1 November 2017 and 1 November 2018, collecting all SNPP CrIS spectra measured on days with maximum VCD less than 1 DU SO₂ present anywhere in the atmosphere (Fig. 3a).

Figure 4 SNPP CrIS mean (a) and standard deviation (b) brightness temperature at 1300 cm⁻¹ for each grid cell and each season interval. The red dot corresponds to the data shown in panel (c). (c) Marginal PDF of the background spectrum indicated by the red dot in panels (a) and (b) with mean spectrum (red dashes) and several individual marginal PDFs (right) shown.

This leaves a database of more than 3.6×10⁸ SO₂-free CrIS spectra over the 1-year period. We classify the spectra regionally and seasonally, partitioning this database into four seasons and $\mathrm{5}{}^{\circ }\times \mathrm{5}{}^{\circ }$ latitude and longitude grid cells yielding 10 368 bins (Fig. 4). Each bin has a full set of marginal PDFs (brightness temperature histograms) for each channel and a covariance matrix characterizing the correlation structure among the channels. This partitioning reduces the overall variability represented in the mean spectra and thus also reduces the magnitude of the error covariance matrix entries while still capturing the fundamental variability due to spectral trends between regions and throughout a year.

For each season–latitude–longitude bin, we construct a representative sample of 10 000 possible background spectra that conform to the set of channel-wise marginal distributions and the covariance matrix relevant to that bin. Although our database is large enough to construct this sample for each bin, we generate these possible spectra via another method because it is preferable (from a mathematical standpoint) that the samples represent only what is known statistically about the spectra. That is a subtle point, but because the channel marginal distributions are represented as histograms (with finite range), generating synthetic background spectra very slightly damps the possible variance contained in a set of real measured spectra and limits the possibility of two key issues: (i) that real but anomalous or erroneous background spectra will be used and (ii) that real spectra with SO₂ just below 1 DU (the database threshold) will be used as a supposedly SO₂-free background. For example, if extreme record-breaking conditions (e.g., hurricanes, droughts, etc.) appear in the database collection interval, their spectra will get caught in the database. These events will affect the covariance matrix and marginal distributions; however, they will cause far more variance in the retrieval if used as backgrounds than if they can only affect the used backgrounds as a forcing on the bin statistics. Additionally, construction of similar databases for other sensors could still proceed with fewer database entries since the statistics of the season–latitude–longitude bins would be expected to converge after fewer entries than were used here.

Since the channel-wise marginal distributions are generally non-Gaussian (e.g., Fig. 4), sampling the random vector Y_bg with covariance matrix S is non-trivial. The general problem of sampling a correlated random vector with known non-normal marginal distributions and covariance matrix is accomplished by a transform sampling technique known as NORTA (NORmal To Anything) (Cario and Nelson, 1997). The NORTA process by which we generate samples of Y_bg is detailed in Appendix A. In the above retrieval (in the pre-screening and fully probabilistic phases), the background spectrum and covariance are interpolated spatially from the collection of binned backgrounds to each CrIS FOV center by a bilinear interpolation scheme using the four nearest season–latitude–longitude cells (Appendix C).

NOAA-20 CrIS has very similar radiometric characteristics as SNPP, except that NOAA-20 FOV 7 noise is within specification and is therefore considered in this study (JPSS CrIS SDR Team, https://www.star.nesdis.noaa.gov/jpss/documents/AMM/N20/CrIS_SDR_Validated.pdf, last access: 24 August 2020). Consequently, we use the SNPP-generated backgrounds in SO₂ retrievals with both SNPP and NOAA-20 CrIS spectra.

3 Results

3.1 Test case I: Raikoke, Kuril Islands, 2019

At approximately 18:00 UTC on 21 June 2019 (04:00 LT), Raikoke volcano in the Kuril Islands erupted for the first time since 1924 (Sennert, 2019; Global Volcanism Program; Hedelt et al., 2019). The strongest pulses of the eruption rose to an altitude of approximately 13 km, forming an umbrella cloud that was quickly advected to the east by strong winds. In the first hours of the eruption, SO₂ columns with VCD>900 DU (Hedelt et al., 2019, >1000 DU, Simon Carn, personal communication, 2019) were detected. The strongest individual measurement made by our method (48.52107^∘ N, 167.25615^∘ W; 15:22:25 UTC, 22 June 2019) had a mean total VCD of 432 DU with a standard deviation of 15 DU (Fig. 5); however, because there is significantly greater uncertainty within the support of the height PDF, the largest value of mean plus uncertainty occurs just below the upper end of the height PDF support (Fig. 5c). This underestimate early in the Raikoke cloud history was likely due two factors, channel saturation despite the specialized strong column retrieval and the fact that the footprint and layout of CrIS FOV leaves many gaps (≈30 % by area) where extremal values could have been present. Because this analysis focuses on the SO₂ ν₃ band (1300–1410 cm⁻¹), it is very unlikely that ash was responsible for the strong underestimation (e.g., Carboni et al., 2012). As described above, more specialized retrievals could be devised to increase the sensitivity to very strong SO₂ loading as in the Raikoke case; however, such schemes are beyond the scope of this work, which is principally concerned with advances in height information and estimating VCD values in more typical (low to moderate concentration) emissions. Within about 1 d, the ash and SO₂ were entrained into a large extratropical cyclone, which heavily distorted the dispersion of the cloud, with the SO₂ cloud being pushed to the north and dispersing in both easterly and westerly directions (Fig. 6). Early in this complex dispersion, SO₂ VCD values remained strong despite a rapid decline in eruptive output. This is most likely a result of the convergence caused by entrainment into the cyclone. Based on the probabilistic retrieval and tropopause data from the National Centers for Environmental Prediction (NCEP), it is clear that the vast majority of the SO₂ cloud mass was in the lower stratosphere, with only a small lower layer in the mid-upper troposphere which had mostly dispersed after the first week (Fig. 6b–l). After 1 month, the SO₂ cloud had spread out over most of the Northern Hemisphere above 30^∘ N with most VCD values <2 DU; however, some columns remained as strong as 20 DU. After 2 months, traces of the SO₂ cloud remained over northern Canada and Hudson Bay with all measured VCD less than 1 DU.

Figure 5 NOAA-20 CrIS mean total VCD (a) and median height (b) early in the evolution of the Raikoke eruption cloud. The star indicates the location of panel (c). (c) Probabilistic retrieval of SO₂ layer altitude (PDF, green) and partial VCD (mean, solid blue; mean ± standard deviation, dotted blue; PDF, color bar) for the strongest individual VCD measured by CrIS in the Raikoke cloud (48.52107^∘ N, 167.25615^∘ W; 22 June 2019, 15:22:25 UTC).

Figure 6 Time evolution (top to bottom) of the Raikoke SO₂ plume from NOAA-20 CrIS showing the expected-value total VCD (a, e, i), the expected-value stratospheric partial VCD (b, f, j), the probability that the SO₂ layer is in the stratosphere (c, g, k), and the median layer height (d, h, l). The height of the tropopause was calculated from daily NCEP reanalysis data.

3.2 Test case II: early detection of SO₂ emission from Bogoslof Volcano, Aleutian Islands, 2016

In the 2016–2017 eruptive period at Bogoslof volcano, 70 explosive events were identified (Coombs et al., 2018, 2019). The first five explosions were not detected in real time and could only be identified and characterized after reanalysis of satellite and other data sources (Coombs et al., 2019). The first CrIS detection of the SO₂ cloud from this sequence of explosions occurred at 22:48 UTC on 16 December 2016 (cluster of 17 CrIS FOVs, approximately 300 km NE of Bogoslof), which was most likely the Event 4 (18:39 UTC) SO₂ plume drifting downwind (Coombs et al., 2018, 2019). As noted in Coombs et al. (2018), the USGS Alaska Volcano Observatory (AVO) was not able to issue a Volcanic Activity Notice (VAN) for this small event and consequently no height information was generated until the reanalyses of Schneider et al. (2020), in which a cloud height of 6.1 km was determined. The SACS near-real-time retrieval (https://sacs.aeronomie.be/) only detected SO₂ from this cloud in two IASI FOVs, which was not sufficient to trigger an alert notification. This small pulse was observable by neither the multispectral infrared remote sensing methods nor automated analysis of multispectral signatures and cloud growth rates (Pavolonis et al., 2013, 2015a, b, 2018; Schneider et al., 2020). CrIS median heights are mainly clustered between 5–8 km with some scatter due to localized cloud edge effects (Fig. 7c, d). This is broadly consistent with the reanalysis of Schneider et al. (2020).

Figure 7 Initial (classical) z score (a), mean total VCD (b), and median height (c) for an explosion from Bogoslof very early in the 2016–2017 eruption (SNPP CrIS). Height PDFs (d) and expected (mean) cloud concentration (e) profiles for the detected SO₂ cloud. Retrievals from the high-noise SNPP CrIS FOV 7 (dotted red lines in d and e) are not shown in panels (a)–(c).

Of particular importance in this small cloud made up of only a few (17) FOVs, SNPP CrIS FOV 7 is significantly nosier (above specification) than that of other FOVs (Zavyalov et al., 2013; Han et al., 2013), and consequently the FOV 7 retrievals are highly suspect and have not been used to estimate the SO₂ height. They are included in Fig. 7d and e mainly to illustrate how an increase in instrument noise affects height information. As might be expected, increased instrument noise propagates larger uncertainty to the height PDFs and tends to distribute their centers to the lower and upper end of the range of considered altitudes. This is the signature of the central role the arg max operation plays in the retrieval. For example, if higher noise is propagated to the z-score height functions, the arg max can produce wildly different heights even for small differences in the z-score height profile.

Because the probabilistic framework allows the calculation of a mean partial VCD, we may derive a formula for the mean or expected concentration profile by similar means as those for Eq. (15):

$$\begin{array}{}\text{(18)}& \mathbb{E}\left[C\right(h\left)\right]=\mathbb{E}\left[{\displaystyle \frac{\mathrm{d}}{\mathrm{d}h}}\widehat{X}\left(h\right)\right]=f_H\left(h\right)\mathbb{E}\left[\widehat{X}|H=h\right],\end{array}$$

which is shown for FOVs in the detected Bogoslof cloud (Fig. 7e). This example demonstrates that the CrIS SO₂ detection and characterization scheme is sufficiently sensitive to capture some small emissions that are generally difficult to observe by other means.

3.3 Test case III: resolving strong stratification in an SO₂ plume, Sierra Negra, Galapagos Islands, 2018

On 26 June 2018, after a period of elevated seismicity, the onset of a major eruption at Sierra Negra was signaled by volcanic tremor at 19:40 UTC, producing an ash and SO₂ plume at 20:09 UTC (Carn et al., 2018; Vasconez et al., 2018; Hedelt et al., 2019). The first CrIS observation also occurred at 20:09 UTC from SNPP, detecting SO₂ above the Sierra Negra on three adjacent FOVs on the edge of scan with maximum initial z scores of approximately 9, 16, and 31. Subsequent overpasses show the plume rising to approximately 14–19 km and spreading in a complex manner due to vertical wind shear, as evidenced by the lower plume altitudes spreading towards the west and the upper plume altitudes spreading towards the southeast (Fig. 8). The significant shearing of the eruption plume enables direct observations of the cloud at many levels. Consequently, this eruption forms a good opportunity to highlight the broad sensitivity of this method in detecting and characterizing SO₂ at every elevation from the vent (1.124 m a.s.l.) up to ∼14–19 km and potentially higher. Additionally, this example highlights the strength of the probabilistic height retrieval, enabling the retrieval of any desired confidence interval on the height. Here we compute the 90 % confidence interval (Fig. 8), highlighting the fact that the 95th and 5th percentiles are in general not symmetric around the median or the same size at different locations within the plume. Because our method retrieves consistent statistics across all measurements, we can ensure the stability of the method and derived probabilities in particular, giving good smoothness even without post-processing. As described above, this is not necessarily the case for other height retrievals, which compute a single estimate with constant uncertainty.

Figure 8 Time evolution (top to bottom) of the 27 June 2018 Sierra Negra SO₂ plume height represented as the 5th percentile (a, d, g), median (b, e, h), and 95th percentile (c, f, i) heights. Data are merged SNPP CrIS and NOAA-20 CrIS with SNPP CrIS FOV 7 excluded.

4 Discussion

4.1 Comparison with other data

Although a deep analysis of the differences between the present method and others is beyond the scope of the present work, here we highlight a brief representative comparison of our SO₂ retrievals with data from TROPOMI, IASI, and the Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP) during the evolution of the Raikoke eruption cloud. These sources of data are described briefly here.

Figure 9 Representative comparison between TROPOMI, IASI, and CrIS SO₂ data. TROPOMI VCD given a 1 km thick layer at 1 km (a), 7 km (b), and 15 km (c) altitudes; IASI (METOP-B) VCD given a 1 km thick layer at the IASI height estimate (f); and CrIS (NOAA-20) mean total VCD (integrated against the height PDF). IASI height estimate (f); CrIS 5th (g), 50th (median, h), and 95th (i) percentile heights; and CrIS initial z score (j). The CrIS–TROPOMI comparison regions analyzed in Fig. 10 are shown in panels (c) and (e) as outlined dashed red and blue latitude–longitude boxes.

As mentioned above, CrIS is a very similar instrument to IASI (both Fourier transform Michelson interferometers). The relevant instrument differences are that IASI (aboard EUMETSAT satellites METOP-A/B, 21:30 LTAN) covers both the ν₁ and ν₃ SO₂ absorption bands (though only ν₃ is used to generate IASI heights), IASI's spectral resolution is 0.5 cm⁻¹ (apodized spectra) compared with 0.625 cm⁻¹ for apodized FSR CrIS, and, as mentioned above, despite IASI's slightly smaller FOV size (12 km-diameter at nadir) compared with CrIS (14 km-diameter at nadir), the greater spatial number density of CrIS FOVs gives CrIS a higher resolution than IASI by virtue of a smaller average sampling distance. For comparison with the set of NOAA-20 CrIS overpasses of the Raikoke cloud in Fig. 9 (24 June 2019, 22:42 UTC, to 25 June 2019, 02:36 UTC), the nearest IASI height data were collected from IASI instrument aboard METOP-B between 00:00 and 11:59 UTC on 25 June 2019. Although some of the data in this interval are highly asynchronous, the cloud did not experience major changes in this time, and the heights (and to a lesser extent VCD) can be compared. As described in Table 1, IASI heights are generated by the Clarisse et al. (2014) retrieval and the VCDs are computed by a related nonlinear technique using the retrieved heights (Clarisse et al., 2012). In the framework presented here, this VCD is the conditional VCD sampled at the IASI-retrieved height.

Figure 10 Comparison between TROPOMI and CrIS SO₂ data for the red (a–e) and blue (f–j) regions outlined in Fig. 9. (a) Relative difference between CrIS FOVs and TROPOMI-synthesized CrIS FOVs for the red-outlined region. (b) The same as panel (a) but smoothed onto a 100 km × 100 km grid. (c) Direct (dots) and rank-order (lines) comparison of CrIS and TROPOMI VCD from (a; grey) and (b; red). (d) Log scale version of panel (c). (e) Histogram (relative frequency) comparison for CrIS and TROPOMI data in red-outlined region, with dotted lines corresponding to the relative frequency of a single measurement. (f–g) The same as panels (a)–(e) for the blue-outlined region.

TROPOMI is a UV spectrometer operating aboard the TROPOMI the Copernicus Sentinel-5 Precursor (S5P) satellite, which orbits only 3.5 min behind SNPP CrIS (Veefkind et al., 2012). TROPOMI represents a significant advance in monitoring of SO₂ and other atmospheric constituents due to its very high spatial resolution (7×3.5 km² pixels at nadir) increasing the sensitivity and signal to noise ratio in a given region by a factor of 3 over OMI and OPMS (Theys et al., 2017, 2019; Fioletov et al., 2020). The TROPOMI SO₂ data used here are generated from back-scattered UV radiances with the S5P operational processing algorithm (a differential optical absorption spectroscopy, DOAS, technique)(Theys et al., 2017). Because most UV-based techniques are not able to retrieve height information, UV-based VCD data are typically presented as conditional VCD for several heights. The TROPOMI SO₂ conditional VCDs are given for assumed SO₂ plume altitudes of 1, 7, and 15 km. As the vast majority of the Raikoke cloud was in the stratosphere, we use the 15 km data for direct comparison with CrIS (Fig. 10). Obviously such high-resolution data are difficult to compare directly with CrIS, so it was first resampled to the CrIS FOV footprints by constructing weighted averages of TROPOMI pixel data with weights determined by the fraction of intersecting area between the pixels and each elliptical CrIS FOV, as in Sun et al. (2018). Unfortunately, since SNPP CrIS was experiencing a major outage during the Raikoke cloud's evolution due to an electrical fault, only NOAA-20 CrIS was able to measure the cloud, and consequently the NOAA-20 CrIS-TROPOMI comparison is asynchronous by about 50–55 min instead of the 3.5–5 min that would have been possible if SNPP CrIS was functioning.

Figure 11 Representative example comparison of CALIOP lidar backscatter (a, b), CrIS SO₂ height PDFs (b partially transparent, c), and CALIOP vertical feature mask (d) for the Raikoke cloud on 25 June 2019. (e) Nearest-neighbor gridded interpolation of CrIS SO₂ median heights with the closest CALIOP overpass shown (black arrow indicating descending orbit <15 min after CrIS acquisition).

The last source of comparison data is 532 nm backscattered lidar measurements from the CALIOP overpass of the Raikoke cloud between 14:32 and 14:36 UTC on 25 June 2019 (Fig. 11). CALIOP profiles aerosols, clouds, and other features between the ground and lower stratosphere (Winker et al., 2009). Although it cannot detect molecular SO₂, it can detect volcanic ash and sulfate aerosols, which are a photochemical product of SO₂ in the atmosphere. As is typical of volcanic SO₂ studies, highly attenuating CALIOP aerosol layers (especially in the stratosphere) are considered here as a proxy for the presence of SO₂ (Clarisse et al., 2014; Carboni et al., 2016, e.g.,).

As mentioned above, our strongest total VCD measurement from the Raikoke cloud was 432 DU. This is significantly lower than the maximum detected by TROPOMI (>900 DU, Hedelt et al., 2019) and several other UV-based methods (∼1000 DU, Simon Carn, personal communication, 2019). This suggests that our method, despite the integrated height estimate and the specialized retrieval for strong SO₂ loading, currently cannot fully capture these extremely high VCD values; however, away from these extreme values, our retrieval performs well in comparison to TROPOMI and IASI (Figs. 9a–e, 10a–e). Other than the relatively few columns with unusually large VCD values, the largest discrepancy between the TROPOMI-retrieved Raikoke cloud and that from CrIS is that the CrIS retrieval does not fully resolve the long, narrow, and diffuse east–west-running cloud to the south of the main cloud, although the CrIS retrieval does resolve some similar features elsewhere in the cloud. This may be due to the very high spatial resolution and sensitivity of the TROPOMI data compared with the coarse CrIS resolution that also contains ≈30 % gaps. An additional contributing factor is that the CrIS retrieval starts with an initial detection of the z score for pre-screening (Fig. 9j) and only retrieves the height and VCD for FOVs with z>5. Traces of this narrow, diffuse cloud are present in the initial z-score field, although it mostly presents with a z score below the detection threshold but slightly above the background noise. Because the western part of this cloud is low altitude, the weak signal in CrIS and IASI but good detection by TROPOMI is evidence that the IR sensors were limited somewhat by interference with larger quantities of water vapor. Overall, it is likely that all of these factors played a role.

A more detailed view of the comparison between CrIS VCDs and TROPOMI conditional VCDs at 15 km altitude is shown for two regions of the Raikoke cloud in Fig. 10. As described above, these measurements are asynchronous by about 50–55 min, which is evident in Fig. 10a, b, f, g, where the largest errors between CrIS and the CrIS-resolution TROPOMI occur at the cloud edges (as well as internally in areas with more complex motion) where it had clearly moved between these observations. These two regions were selected because they both contain VCDs spanning several orders of magnitude but different ranges. The northwestern region (dashed red, Figs. 9, 10a–e) contains VCD values generally up to about 50 DU whereas the southeastern region (dashed blue, Figs. 9, 10f–j) contains VCD values as high as 230 DU. As can be seen from the histograms of native resolution and CrIS-resolution TROPOMI, the lower spatial resolution (and gaps) of CrIS is at least partially responsible for its inability to resolve some of the highest VCD values measured by TROPOMI in native resolution. In regions of the cloud with generally lower concentrations, CrIS and TROPOMI compare well, scaling approximately one-to-one, with much of the noise being generated by the asynchrony. Although not a perfect comparison, the rank-order correlation is also shown for the CrIS-resolution comparison and a coarse, 100 km pixel aggregate. Taken together, the FOV-by-FOV and 100 km pixel comparisons (both exact matching and rank-order) demonstrate that CrIS matches TROPOMI VCDs moderately well for these low to moderate concentration clouds. In higher concentration clouds (Fig. 10f–j), the comparison has a similar noise profile; however, CrIS consistently underestimates TROPOMI by about 75 % (−0.25 relative error, Fig. 10f, g). Based on the histograms and correlation of CrIS and CrIS-resolution TROPOMI, there is a similar distribution of VCDs below about 50 DU that becomes divergent by about 100 DU. In all regions of the Raikoke cloud, the sensitivity of the fully probabilistic retrieval is limited to VCDs greater than about 0.3–0.5 DU.

The main focus and strength of our approach is the ability to generate physics-based PDFs for the height. Because our retrieval is based on the operational algorithm in use for IASI, our retrieved heights are very similar to those from IASI, although there are key differences readily apparent in Fig. 9f–i. As mentioned above, the IASI heights represent the height retrieved due to the mean background spectrum; however, because the retrieval of height is not linear (due to the arg max operation), the retrieved height is not the expected value height. Inspection of the PDFs generated by this approach show that they are typically non-symmetric, although exact comparison is not possible due the orbital separation between the satellites carrying IASI (METOP-A,B) and those carrying CrIS (SNPP, NOAA-20). At least for the Raikoke cloud, the IASI heights are almost entirely bound within the CrIS 90 % confidence interval (Fig. 9f, g, i) and are similar to the CrIS median heights (Fig. 9h) in most areas. For the snapshot of the Raikoke cloud shown in Fig. 9, the largest differences in height appear at the southernmost part of the cloud, where the IASI heights are >19 km in altitude over a significant region. Furthermore, the IASI height estimate (Fig. 9f) varies significantly over nearby continuous parts of the cloud, whereas the CrIS median height is more consistent across space with some minor variation due to noise.

Although not shown here, Hedelt et al. (2019) have recently developed a new SO₂ height retrieval for TROPOMI using inverse learning machines. Although computationally expensive to train, such an approach has the advantage of computation speed of the inversion once deployed, though it has not yet been incorporated into the TROPOMI SO₂ data product as of this writing, and thus a direct comparison is not possible. However, it is clear from the data presented in Hedelt et al. (2019) that their height estimates for the Raikoke cloud span the CrIS 90 % confidence interval, though they are significantly nosier than CrIS SO₂ height estimates and display a very prominent negative trend in height versus VCD, leading to systematically higher layer heights on the cloud edges than in the cloud centers (Fig. 14 of Hedelt et al., 2019). CrIS PDFs very rarely show a similar trend.

Because we retrieve a PDF on each CrIS FOV rather than a single estimate, we can compare the PDFs directly to data from a CALIOP overpass of the cloud. Here we show an example comparison from Raikoke; however, a full comparison for every overpass of the Raikoke cloud is the subject of future work. For the first several days after the eruption, there was still significant ash suspended in the dispersing cloud, leading to the appearance of several highly attenuating layers in CALIOP data between 10–15 km (Fig. 11a, b). The comparison we focus on is between CrIS retrievals from 14:18:00 to 14:24:00 UTC and CALIOP data from a subsequent overpass between 14:32:06 and 14:36:14 UTC on 25 June 2019. To prevent the mixing of nearby height PDFs, the CrIS data are interpolated to the CALIOP track by nearest-neighbor interpolation, creating a profile of the nearest CrIS SO₂ height PDFs. The nearest-neighbor interpolation over all of space is shown in Fig. 11e with 16 km pixels (the average CrIS sampling distance).

Overall, there is good agreement between the CrIS SO₂ height PDFs and the altitudes of the strongly attenuating CALIOP layers; however, the CrIS PDFs have several important characteristics that complicate comparison. There is some minor noise derived mainly from two anomalous CrIS FOVs at the cloud edge ($\approx \mathrm{62}{}^{\circ }$ N, 175^∘ W) exactly at the CALIOP track, leading to anomalously high altitudes there (Fig. 11e). Most interestingly, some regions of the cloud have bimodal PDFs (Fig. 11b, c, south of 60^∘ N). Preliminary investigations of this retrieval suggest that such occurrences are not exceedingly rare throughout this and other clouds; however, such PDFs occur widely and tend to persist over moderate distances (as in Fig. 11b, c).

In the strictest sense, such PDFs can only be attributed to the presence of similar statistical features in the background. Specifically, if the background spectrum probability space is dominated by two sets of meteorological conditions (for example, one mode representing deep convective cloud radiances and another for cloud-free radiances), then multiple populations of the Monte Carlo height samples may accumulate, leading to a multimodal height PDF.

Relaxing this strict interpretation, there is some evidence that these bimodal PDFs may represent the presence of SO₂ at multiple altitudes. In this case, the strongest CrIS probabilities occur in a lower layer around 12 km altitude (Fig. 11b, c) suggesting the presence of molecular SO₂ layer; however, there is no attenuating CALIOP layer there. Additionally, the CrIS retrieval assigns significant (but less) probability mass to a higher level colocated with a strongly attenuating CALIOP layer at about 15 km. Since these data are very early in the cloud's evolution (<5 d), this layer is likely dominated by volcanic ash rather than being dominated by sulfate aerosols, having not had enough time to convert large portions of the erupted SO₂. Of note, because the upper CALIOP layer is very strongly attenuating, it may completely shadow any evidence of lower diffuse ash clouds if they existed. Considering that ash minimally affects the SO₂ ν₃ band, this colocation at 15 km altitude suggests that this strongly attenuating layer (likely mainly ash) contains SO₂. The fact that the stronger of the two probabilities is in the lower layer combined with the CALIOP–CrIS colocation could be interpreted as evidence that there is SO₂ at both layers. Although this is a possibility, confirmation of such a configuration would require a deeper analysis including forward and inverse trajectory modeling with advanced data assimilation techniques and therefore is well beyond the scope of the present work.

4.2 Long-term analysis of the Raikoke SO₂ cloud

By retrieving PDFs for height and partial VCD it is possible to enhance time series analysis of SO₂ clouds accordingly, enabling the generation of time series with quantified uncertainty. As an example, we calculate the total and stratospheric SO₂ mass time series probabilistically as sums of many independent retrievals (Fig. 12).

Figure 12 (a) Total (blue) and stratospheric (red) mass of SO₂ above 30^∘ N for 2 months following the eruption of Raikoke showing CrIS-derived mean (solid lines) ±1 standard deviation (dotted lines), with the log scale in the inset. Black shows the preliminary TROPOMI mass (Fig. 12; https://so2.gsfc.nasa.gov/pix/special/2019/raikoke/raikoke_tropomi_so2stl_0621-07062019.html, NASA/GSFC, 2019). “Instantaneous” (daily aggregated) e-folding time PDFs for the total (b) and stratospheric (c) SO₂ in the Raikoke eruption cloud.

To estimate the mass of an SO₂ cloud, the values of the retrieval on the CrIS FOVs must be interpolated onto a continuous grid spanning the cloud. In the present study we use an equal-area grid and perform nearest-neighbor interpolation of the CrIS FOV retrievals. We calculate the mass (with uncertainty) as in Appendix E for both the total SO₂ mass in the atmosphere and in the stratosphere only. The stratospheric partial VCD values $\widehat{X}\left(h_{\mathrm{tropopause}}\right)$ were calculated from Eqs. (17a) and (17b) using daily NCEP/NCAR tropopause pressure level and geopotential height reanalysis data (Kalnay et al., 1996). By exploiting the fact that the mass is a normal random variable at each time, we can from these results estimate the daily “instantaneous” apparent e-folding time of the SO₂ as $\mathit{\tau }=-M\left(t\right)/\dot{M}\left(t\right)$ with $\dot{M}\left(t\right)$ calculated by finite difference and the PDFs computed by standard results in probability theory (Fieller, 1932; Hinkley, 1969; DeGroot and Schervish, 2012, Appendix E). The e-folding times shown here are really apparent measurements since mass is lost from the cloud not only due to photochemical conversion of SO₂ to sulfate aerosols but also due to the dispersion and dilution of the cloud below levels that can be detected. We include this time series because it illustrates some interesting aspects of the Raikoke cloud's evolution and because it illustrates a logical extension of how uncertainty is propagated through this work.

From Fig. 12 it is clear that the SO₂ cloud, as characterized by CrIS, did not immediately show the exponential mass decay that has been used in similar studies of large eruptions (e.g., Read et al., 1993; Carn et al., 2017; Krotkov et al., 2010). The observed trend is likely a combination of retrieval limitations and genuine atmospheric chemical properties. As described above, the CrIS VCD underestimate early in the Raikoke cloud history was likely due to channel saturation despite the specialized strong column retrieval. The strongest CrIS VCDs were only approximately 50 % of the strongest columns reported at the time; 432 DU from CrIS versus >900 DU (Hedelt et al., 2019) and >1000 DU (Simon Carn, personal communication, 2019). In aggregate, CrIS total VCD were approximately 75 % of those of TROPOMI in the most concentrated regions (Fig. 10). This underestimate was propagated to the CrIS maximum total mass estimate of approximately 1.1 Tg SO₂ (Fig. 12) compared with preliminary estimates from other sources (1.4–1.5 Tg; Sennert, 2019; Global Volcanism Program; Simon Carn, personal communication, 2019). Very large VCD values persisted for many days, though after approximately 10 d most of the cloud had diluted sufficiently that CrIS and preliminary TROPOMI data were in good agreement between 10 and 20 d from the eruption, both decreasing from about 1.0 to 0.6 Tg (Fig. 12; https://so2.gsfc.nasa.gov/pix/special/2019/raikoke/raikoke_tropomi_so2stl_0621-07062019.html, NASA/GSFC, 2019).

Despite this early underestimation and the fact that the SO₂ injection is more or less instantaneous compared with the time span of cloud detectability, CrIS total mass (and, to a lesser extent, preliminary TROPOMI total mass as well) did not begin the anticipated exponential decay until approximately 15–20 d after SO₂ injection. We posit that this was mainly due to CrIS channel saturation and highly irregular dispersal by extratropical cyclonic winds affecting the total CrIS-measurable mass and cloud dilution, respectively. However, the presence of this effect in preliminary TROPOMI suggests that at least some contribution may have derived from time-dependent SO₂ to sulfate conversion kinetics. As has been invoked in previous studies, this could have been the result of limited hydroxyl radical (OH) early in the cloud history (e.g., Theys et al., 2013; Sekiya et al., 2016), although detailed chemical modeling to support this is beyond the scope of this work. As the e-folding times shown in Fig. 12c and d are calculated by finite difference, they are quite noisy; however, they exhibit the same trend of slow decay early in the cloud history even after CrIS and TROPOMI masses become similar. After about 35 d, the apparent e-folding time for the total VCD settles to a more or less constant value (∼10 d e-folding time). The apparent stratospheric e-folding time exhibits a very different pattern, remaining approximately constant (∼10 d e-folding time) until about day 30, when it begins to shorten. The large uncertainty in the stratospheric e-folding time after approximately 30 d is the result of the increasing relative stratospheric mass uncertainty and the particular details of the PDF calculation (Hinkley, 1969). Because SO₂ is typically assumed to have a long lifetime in the stratosphere, this is likely the result of dilution below the detection threshold over large low-concentration regions of the cloud, although some portion of the SO₂ re-entering the troposphere as it extends to lower latitudes (where the tropopause is higher) may also play a role considering that the total column does not show this effect as significantly.

5 Conclusions

• i.

  New probabilistic enhancement of existing hyperspectral IR SO₂ retrieval techniques enables the retrieval of PDFs for SO₂ height and partial VCD, providing significant statistical power, precision, and consistency to global SO₂ detection, tracking, and characterization efforts. Retrieving these PDFs enables the calculation of many new quantities, including exceedance probabilities for concentration, layer height constraint probabilities, mean concentration profiles, and mass at different layers with uncertainty. Although these capabilities are primarily beneficial for operational SO₂ monitoring, these methods are relevant to climatological studies because of the ability to the stratospheric fraction of total mass for a given SO₂ cloud if tropopause heights are available from ancillary sources.

• ii.

  This technique is capable of resolving larger VCD values than would be anticipated for a linearized approach due to two factors: (i) the use of the height PDF increases the retrieved total VCD compared with a single height estimate and (ii) the use of a specialized channel subset retrieval that improves the linear approximation when the signal is certain to be dominated by SO₂. However, these techniques are still limited in their ability to resolve extreme VCD values (many hundreds of Dobson units) like some of those observed in the recent (2019) eruption of Raikoke Volcano, Kuril Islands. Because of the improved spatial resolution over IASI and the technique's sensitivity, we can resolve heights for small clouds that cannot be resolved well by IASI. Additionally, the technique can adequately resolve height information across a broad range of plume altitudes, including those in the lower stratosphere and close to the surface (though with more limited detection).

• iii.

  Preliminary comparisons suggest that this method generally compares well with other measurements of SO₂ VCD and altitude; however, the probabilistic framework adds significant value, especially in the retrieval of height information. Cross sections through these probability clouds compare very well with cloud heights from CALIOP lidar backscatter.

• iv.

  As a logical extension of the probabilistic framework, this technique enables the characterization of SO₂ clouds for long-term probabilistic analyses of cloud evolution and key time-varying parameters such as total mass, stratospheric mass, and apparent e-folding time.

• v.

  The algorithms presented here are currently being integrated into VOLCAT, where they will be used for operational SO₂ cloud detection, characterization, and tracking in support of aviation safety. We anticipate future work to include more comprehensive comparison of height PDFs with CALIOP lidar backscatter data, application of these techniques to similar instruments including IASI and the Atmospheric Infrared Sounder (AIRS), and the development of volcanic degassing and aviation-focused products.

Appendix A: Generating correlated random spectra for Monte Carlo retrieval: NORTA sampling

The general problem of sampling a correlated random vector ($Y_i\in \mathit{Y}\in {\mathbb{R}}^N$) with non-normal marginal distributions $F_{Y_i}$ and covariance matrix S is accomplished by a transform sampling technique known as NORTA (NORmal To Anything) (Cario and Nelson, 1997), in which the desired random vector is written as a component-wise inverse marginal transform of a standard normal random vector:

$$\begin{array}{}\text{(A1)}& Y_i=F_{Y_i}^{-\mathrm{1}}\left(\mathrm{\Phi }\left(Z_i\right)\right),\end{array}$$

where Φ is the univariate standard normal cumulative density function (CDF). In this formulation, the individual marginals can be used to transform the components of the standard normal vector. The crux of this method lies in generating a standard normal random vector Z with appropriate correlation structure ρ_Z, which upon component-wise transformation, as above, will produce the desired random vector Y with the desired covariance structure S. That is, a correlated Z must be generated with the correlation matrix

$$\begin{array}{}\text{(A2)}& {\mathit{\rho }}_{\mathit{Z}}=\mathrm{Cov}\left(\mathit{Z}\right)=\mathbb{E}\left[\mathit{Z}{\mathit{Z}}^{\mathrm{T}}\right],\end{array}$$

which, after transformation, results in the correlation matrix ρ_Y=Corr(Y) associated with the known covariance matrix S.

This is accomplished by solving

$$\begin{array}{}\text{(A3)}& {\mathit{\rho }}_{\mathit{Y}}(i,j)=C_{ij}\left[{\mathit{\rho }}_{\mathit{Z}}(i,j)\right],\end{array}$$

for each of the unique $N(N-\mathrm{1})/\mathrm{2}$ correlations (ρ_Z(i,j)) in the lower triangle of ρ_Z. The correlation transformation function is

$$\begin{array}{}\text{(A4)}& C_{ij}\left[{\mathit{\rho }}_{\mathit{Z}}(i,j)\right]:={\displaystyle \frac{\mathbb{E}\left[Y_i\right(Z_i\left)Y_j\right(Z_j\left)\right]-\mathbb{E}\left[Y_i\right]\mathbb{E}\left[Y_j\right]}{\sqrt{\mathrm{Var}\left[Y_i\right]\mathrm{Var}\left[Y_i\right]}}},\end{array}$$

and

$$\begin{array}{}\text{(A5)}& \begin{array}{rl}\mathbb{E}& \left[Y_i\left(Z_i\right)Y_j\left(Z_j\right)\right]=I\left[{\mathit{\rho }}_{\mathit{Z}}(i,j)\right]:=\\ & \underset{{\mathbb{R}}^{\mathrm{2}}}{\iint }{F}_{Y_i}^{-\mathrm{1}}\left[\mathrm{\Phi }\left(z_i\right)\right]F_{Y_j}^{-\mathrm{1}}\left[\mathrm{\Phi }\left(z_j\right)\right]\mathit{\phi }\left(z_i,z_j;{\mathit{\rho }}_{\mathit{Z}}(i,j)\right)\\ & \mathrm{d}{z}_i\mathrm{d}{z}_j,\end{array}\end{array}$$

where $\mathit{\phi }(z_i,z_j;{\mathit{\rho }}_{\mathit{Z}}(i,j\left)\right)$ is the bivariate standard normal density function with correlation ρ_Z(i,j) between the variables Z_i and Z_j. These problems may be inverted individually by various methods outlined in the original formulation of Cario and Nelson (1997).

A1 Application of NORTA to background spectrum

For the purpose of generating the correlated random background spectrum in the present study, we make several modifications to the classical method of Cario and Nelson (1997) to make our problem tractable but retain the high fidelity of the CrIS measurements. In the present study, 177 CrIS channels are used, representing the FSR mid-wave band between 1300–1410 cm⁻¹. This yields 15 576 independent correlations – matching inverse problems that are solved in the present study by gradient descent iteration. As in the main text, we first collect a database of SO₂-free CrIS spectra for each seasonal $\mathrm{5}{}^{\circ }\times \mathrm{5}{}^{\circ }$ bin and compute the 1300–1410 cm⁻¹ background covariance matrix S, as well as the channel-wise marginal distributions, as a sequence of 177 histograms. Then we perform the NORTA process to generate 10 000 samples of Y_bg for each season–latitude–longitude bin.

A2 Numerical integration of the joint expectation

Because each of the correlations matching inverse problems requires multiple rounds of numerical integration of Eq. (A5), we make several modifications to increase computational efficiency. Because many of the channel correlations are strong, the bivariate standard normal distribution $\mathit{\phi }(z_i,z_j;{\mathit{\rho }}_{\mathit{Z}}(i,j\left)\right)$ for each such pair of channels is very narrow. Consequently, for a typical rectangular sampling domain, the integrand of Eq. (A5) is approximately zero almost everywhere except for a concentrated region in which accurate numerical integration is challenging. To reduce wasted integrand samples, we make a standard transformation of the bivariate normal distribution and then cast the domain in polar coordinates and approximate the integral over a finite radial domain [0,R]:

$$\begin{array}{}\text{(A6)}& I\left[{\mathit{\rho }}_{\mathit{Z}}(i,j)\right]:=\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{R}{\int }}{G}_{ij}\left(r,\mathit{\theta };{\mathit{\rho }}_{\mathit{Z}}(i,j)\right){\displaystyle \frac{e^{-r^{\mathrm{2}}/\mathrm{2}}}{\mathrm{2}\mathit{\pi }}}r\mathrm{d}r\mathrm{d}\mathit{\theta },\end{array}$$

where

$$\begin{array}{}\text{(A7)}& \begin{array}{rl}{G}_{ij}& \left(r,\mathit{\theta };{\mathit{\rho }}_{\mathit{Z}}(i,j)\right)=F_{Y_i}^{-\mathrm{1}}\left[\mathrm{\Phi }\left(r\cos \mathit{\theta }\right)\right]\\ & \times F_{Y_j}^{-\mathrm{1}}\left[\mathrm{\Phi }\left({\mathit{\rho }}_{\mathit{Z}}(i,j)r\cos \mathit{\theta }+\sqrt{\mathrm{1}-{\mathit{\rho }}_{\mathit{Z}}^{\mathrm{2}}(i,j)}r\sin \mathit{\theta }\right)\right],\end{array}\end{array}$$

is the product of the transformed ith and jth components of the desired random vector. This approximates the full improper integral within tolerance ϵ, with the fixed value $\mathit{ϵ}={\mathrm{10}}^{-\mathrm{6}}$ in the present study. The terminal radius is estimated conservatively as

$$\begin{array}{}\text{(A8)}& R=\sqrt{\mathrm{2}\ln \left({\displaystyle \frac{\mathrm{1}}{\mathit{ϵ}}}{\displaystyle \frac{max{Y}_{i}max{Y}_j}{min{Y}_{i}min{Y}_j}}\right)},\end{array}$$

where the maximum and minimum values of the components were recorded during the generation of S and the marginal distributions from sample data. In the absence of knowledge of these minima and maxima, they could be estimated from the marginal distributions as fixed percentiles of these distributions.

The estimated radial limit of integration R is derived by requiring that

$$\begin{array}{}\text{(A9)}& \underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{R}{\overset{\mathrm{\infty }}{\int }}{G}_{ij}{\displaystyle \frac{e^{-r^{\mathrm{2}}/\mathrm{2}}}{\mathrm{2}\mathit{\pi }}}r\mathrm{d}r\mathrm{d}\mathit{\theta }\le \mathit{ϵ}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{G}_{ij}{\displaystyle \frac{e^{-r^{\mathrm{2}}/\mathrm{2}}}{\mathrm{2}\mathit{\pi }}}r\mathrm{d}r\mathrm{d}\mathit{\theta }.\end{array}$$

This inequality can be solved for R under the conservative estimation that $max{G}_{ij}=max{Y}_{i}max{Y}_j$ and $min{G}_{ij}=min{Y}_{i}min{Y}_j$ which are derived from the maximum and minimum brightness temperatures on each channel from among all background spectra in the database. This leads to

$$\begin{array}{}\text{(A10)}& \underset{R}{\overset{\mathrm{\infty }}{\int }}r{e}^{-r^{\mathrm{2}}/\mathrm{2}}\mathrm{d}r\le \mathit{ϵ}{\displaystyle \frac{min{Y}_{i}min{Y}_j}{max{Y}_{i}max{Y}_j}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}r{e}^{-r^{\mathrm{2}}/\mathrm{2}}\mathrm{d}r,\end{array}$$

and, upon integrating, R may be chosen conservatively to satisfy

$$\begin{array}{}\text{(A11)}& e^{-R^{\mathrm{2}}/\mathrm{2}}\le \mathit{ϵ}{\displaystyle \frac{min{Y}_{i}min{Y}_j}{max{Y}_{i}max{Y}_j}}.\end{array}$$

This transformation to a scaled polar coordinates ensures that the curvature and gradients in the integrand are as small as possible, ensuring that numerical integration is as accurate as possible for a given domain sampling up to the order of the method employed. Additionally, the sensible radial limitation of this integral ensures that it can be computed efficiently to within tolerance without the inclusion of samples which minimally affect the total (approximately Gaussian tails).

A3 Numerical solution of inverse problems

Although each channel pair inverse problem can be solved separately by Newton's method or other algorithms, we solve all of the problems jointly, restating the problem as a gradient and descent minimization of the total square correlation error

$$\begin{array}{}\text{(A12)}& {\mathit{\epsilon }}^{\mathrm{2}}\left({\mathit{r}}_{\mathit{Z}}\right)=\left[\mathit{C}\right({\mathit{r}}_{\mathit{Z}})-{\mathit{r}}_{\mathit{Y}}{]}^{\mathrm{T}}[\mathit{C}\left({\mathit{r}}_{\mathit{Z}}\right)-{\mathit{r}}_{\mathit{Y}}],\end{array}$$

where the vectors ${\mathit{r}}_{\mathit{Z}},{\mathit{r}}_{\mathit{Y}}\in {\mathbb{R}}^{\mathrm{15},\mathrm{576}}$ are comprised of the unique lower triangular elements of ρ_Z and ρ_Y, respectively, and $\mathit{C}:{\mathbb{R}}^{\mathrm{15},\mathrm{576}}\to {\mathbb{R}}^{\mathrm{15},\mathrm{576}}$ is the correlation transformation function cast as a vector function. Although this seems intractable due to the extreme dimensionality, there is no correlation between these dimensions (no correlation between pairwise inverse problems), and consequently the Jacobian is zero everywhere except for its main diagonal. The gradient descent method of Barzilai and Borwein (1988) produces fast convergence to a global minimum due to the monotonicity and bounding properties of the correlation functions C_ij as described by Cario and Nelson (1997). At each iteration, the NORTA process is completed and the error in the synthesized channel marginal distributions and spectral covariance matrix are used as convergence criteria. The vectorization of these independent problems allows for standardization in convergence criteria using global (total) error minimization. Because the error is measured on the final product (the desired random spectra), a minimum number of iterations is required by comparison with performing each problem separately and then generating the desired random spectra.

Appendix B: Derivation of VCD mean, variance, and covariance formulae

As in the text, SO₂ is assumed to exist in a narrow, 1 km thick layer, represented by a box profile:

$$\begin{array}{}\text{(B1)}& \mathit{x}=X\mathrm{\Pi }(h-H)=\left\{\begin{array}{ll}X/L& L/\mathrm{2}<h-H<L/\mathrm{2}\\ \mathrm{0}& \text{otherwise}\end{array}\right.,\end{array}$$

where L is 1 km. For the purposes of simplicity in computation, we consider a limiting case of a very thin layer (L→0), where the finite box profile converges in distribution to the Dirac delta:

$$\begin{array}{}\text{(B2)}& \mathrm{\Pi }(h-H)\stackrel{d}{⟶}\mathit{\delta }(h-H)\text{ as }L\to \mathrm{0}.\end{array}$$

The calculation of the mean and variance below makes extensive use of Fubini's theorem, allowing reordering of iterated integration. Because these integrals contain the Dirac delta, it is not simple to show that the conditions of Fubini's theorem are satisfied due to the Lebesgue non-integrability of the Dirac delta, which in turn stems from the fact that the Dirac delta is not a proper function but is a distribution instead. However, we proceed assuming that the iterated integrals can be interchanged. We remark that a proof using limits of functions (“nascent delta functions”) approaching the Dirac delta could be substituted here at great cost to simplicity and readability. For a finite thickness layer (as assumed in the main text), all convolutions with the Dirac delta below would be replaced with convolutions with a boxcar function (a finite nascent delta). For the 1 km thick box profile above, the results of convolution will be smoothed, 1 km running-average versions of the desired functions. The final integral formulas would be the same except that the integrands would first be smoothed by a 1 km running average. In a discrete setting with 1 km sampling of the retrieved height PDF, the assumption of a 1 km thick box profile yields the same results as using the Dirac delta.

The mean partial VCD is calculated as the expectation

$$\begin{array}{}\text{(B3)}& \begin{array}{rl}\mathbb{E}& \left[\widehat{X}\left(h\right)\right]=\mathbb{E}\left[\underset{\mathrm{0}}{\overset{h}{\int }}\widehat{X}\mathit{\delta }(\mathit{\eta }-H)\mathrm{d}\mathit{\eta }\right]\\ & =\underset{\mathcal{H}}{\int }\underset{\mathcal{X}}{\int }\left(\underset{\mathrm{0}}{\overset{h}{\int }}\widehat{x}\mathit{\delta }\left(\mathit{\eta }-h^{\prime }\right)\mathrm{d}\mathit{\eta }\right)f_{\widehat{X},H}\left(\widehat{x},h^{\prime }\right)\mathrm{d}\widehat{x}\mathrm{d}{h}^{\prime }\\ & =\underset{\mathcal{H}}{\int }\underset{\mathcal{X}}{\int }\left[\left(\underset{\mathrm{0}}{\overset{h}{\int }}\widehat{x}\mathit{\delta }\left(\mathit{\eta }-h^{\prime }\right)\mathrm{d}\mathit{\eta }\right)f_{\widehat{X}|H}\left(\widehat{x}|h^{\prime }\right)\right.\\ & \left.f_H\left(h^{\prime }\right)\right]\mathrm{d}\widehat{x}\mathrm{d}{h}^{\prime },\end{array}\end{array}$$

where ℋ and 𝒳 are the domains of height and VCD (both σ finite measure spaces) and the last equality follows from the definition of the conditional PDF. Rearranging this iterated integral gives

$$\begin{array}{}\text{(B4)}& \begin{array}{rl}\mathbb{E}& \left[\widehat{X}\left(h\right)\right]=\underset{\mathrm{0}}{\overset{h}{\int }}\underset{\mathcal{H}}{\int }\left[\mathit{\delta }\left(\mathit{\eta }-h^{\prime }\right)\right.\\ & \left.f_H\left(h^{\prime }\right)\left(\underset{\mathcal{X}}{\int }\widehat{x}{f}_{\widehat{X}|H}\left(\widehat{x}|h^{\prime }\right)\mathrm{d}\widehat{x}\right)\right]\mathrm{d}{h}^{\prime }\mathrm{d}\mathit{\eta }\\ & =\underset{\mathrm{0}}{\overset{h}{\int }}\underset{\mathcal{H}}{\int }\mathit{\delta }\left(\mathit{\eta }-h^{\prime }\right)f_H\left(h^{\prime }\right)\mathbb{E}\left(\widehat{X}|H=h^{\prime }\right)\mathrm{d}{h}^{\prime }\mathrm{d}\mathit{\eta }.\end{array}\end{array}$$

Because of the symmetry of the delta function ($\mathit{\delta }(\mathit{\eta }-h^{\prime })=\mathit{\delta }(h^{\prime }-\mathit{\eta })$), the integral properties of the delta function yields

$$\begin{array}{}\text{(B5)}& \mathbb{E}\left[\widehat{X}\left(h\right)\right]=\underset{\mathrm{0}}{\overset{h}{\int }}{f}_H\left(\mathit{\eta }\right)\mathbb{E}\left(\widehat{X}|H=\mathit{\eta }\right)\mathrm{d}\mathit{\eta }.\end{array}$$

Because the algebraic form of the variance is $\mathrm{Var}\left(\widehat{X}\right(h\left)\right)=\mathbb{E}\left[{\widehat{X}}^{\mathrm{2}}\right(h\left)\right]-\left[\mathbb{E}\right(\widehat{X}\left(h\right)){]}^{\mathrm{2}}$, only the second moment of the partial VCD $\mathbb{E}\left[{\widehat{X}}^{\mathrm{2}}\right(h\left)\right]$ must be calculated to complete the formula. This quantity is

$$\begin{array}{}\text{(B6)}& \begin{array}{rl}\mathbb{E}& \left[{\widehat{X}}^{\mathrm{2}}\left(h\right)\right]=\mathbb{E}\left[{\left(\underset{\mathrm{0}}{\overset{h}{\int }}\widehat{X}\mathit{\delta }(\mathit{\eta }-H)\mathrm{d}\mathit{\eta }\right)}^{\mathrm{2}}\right]\\ & =\underset{\mathcal{H}}{\int }\underset{\mathcal{X}}{\int }\left[\left(\underset{\mathrm{0}}{\overset{h}{\int }}\widehat{x}\mathit{\delta }\left({\mathit{\eta }}_{\mathrm{0}}-h^{\prime }\right)\mathrm{d}{\mathit{\eta }}_{\mathrm{0}}\right)\left(\underset{\mathrm{0}}{\overset{h}{\int }}\widehat{x}\mathit{\delta }\left(\mathit{\eta }-h^{\prime }\right)\mathrm{d}\mathit{\eta }\right)\right.\\ & \left.f_{\widehat{X}|H}\left(\widehat{x}|h^{\prime }\right)f_H\left(h^{\prime }\right)\right]\mathrm{d}\widehat{x}\mathrm{d}{h}^{\prime }\\ & =\underset{\mathrm{0}}{\overset{h}{\int }}\underset{\mathrm{0}}{\overset{h}{\int }}\underset{\mathcal{H}}{\int }\left[\mathit{\delta }\left({\mathit{\eta }}_{\mathrm{0}}-h^{\prime }\right)\mathit{\delta }\left(\mathit{\eta }-h^{\prime }\right)f_H\left(h^{\prime }\right)\right.\\ & \left.\left(\underset{\mathcal{X}}{\int }{\widehat{x}}^{\mathrm{2}}{f}_{\widehat{X}|H}\left(\widehat{x}|h^{\prime }\right)\mathrm{d}\widehat{x}\right)\right]\mathrm{d}{h}^{\prime }\mathrm{d}{\mathit{\eta }}_{\mathrm{0}}\mathrm{d}\mathit{\eta }\\ & =\underset{\mathrm{0}}{\overset{h}{\int }}\underset{\mathrm{0}}{\overset{h}{\int }}\underset{\mathcal{H}}{\int }\left[\mathit{\delta }\left(h^{\prime }-{\mathit{\eta }}_{\mathrm{0}}\right)\mathit{\delta }\left(h^{\prime }-\mathit{\eta }\right)f_H\left(h^{\prime }\right)\right.\\ & \left.\mathbb{E}\left({\widehat{X}}^{\mathrm{2}}|H=h^{\prime }\right)\right]\mathrm{d}{h}^{\prime }\mathrm{d}{\mathit{\eta }}_{\mathrm{0}}\mathrm{d}\mathit{\eta }\\ & =\underset{\mathrm{0}}{\overset{h}{\int }}\underset{\mathrm{0}}{\overset{h}{\int }}\mathit{\delta }\left({\mathit{\eta }}_{\mathrm{0}}-\mathit{\eta }\right)f_H\left({\mathit{\eta }}_{\mathrm{0}}\right)\mathbb{E}\left({\widehat{X}}^{\mathrm{2}}|H={\mathit{\eta }}_{\mathrm{0}}\right)\mathrm{d}{\mathit{\eta }}_{\mathrm{0}}\mathrm{d}\mathit{\eta }.\end{array}\end{array}$$

Since the dummy variable η always runs between $\mathrm{0}<\mathit{\eta }<h$, the delta function δ(η₀−η) is always centered in the interval $\mathrm{0}<{\mathit{\eta }}_{\mathrm{0}}<h$ and the integral properties of the delta function can be applied again:

$$\begin{array}{}\text{(B7)}& \mathbb{E}\left[{\widehat{X}}^{\mathrm{2}}\left(h\right)\right]=\underset{\mathrm{0}}{\overset{h}{\int }}{f}_H\left(\mathit{\eta }\right)\mathbb{E}\left({\widehat{X}}^{\mathrm{2}}|H=\mathit{\eta }\right)\mathrm{d}\mathit{\eta }.\end{array}$$

Substitution of $\mathbb{E}\left({\widehat{X}}^{\mathrm{2}}\right|H=\mathit{\eta })=\mathrm{Var}(\widehat{X}|H=\mathit{\eta })+\left[\mathbb{E}\right(\widehat{X}|H=\mathit{\eta }){]}^{\mathrm{2}}$ into the formula for $\mathbb{E}\left[{\widehat{X}}^{\mathrm{2}}\right(h\left)\right]$ and subsequently into the formula above for $\mathrm{Var}\left(\widehat{X}\right(h\left)\right)$ completes the variance.

Similarly, for the covariance between partial VCD at two altitudes h=a and h=b with b≥a, only the mixed expectation $\mathbb{E}\left[\widehat{X}\right(b\left)\widehat{X}\right(a\left)\right]$ must be calculated:

$$\begin{array}{}\text{(B8)}& \begin{array}{rl}\mathbb{E}& \left[\widehat{X}\left(a\right)\widehat{X}\left(b\right)\right]\\ & =\mathbb{E}\left[\underset{\mathrm{0}}{\overset{a}{\int }}\widehat{X}\mathit{\delta }\left({\mathit{\eta }}_a-H\right)\mathrm{d}{\mathit{\eta }}_a\underset{\mathrm{0}}{\overset{b}{\int }}\widehat{X}\mathit{\delta }\left({\mathit{\eta }}_b-H\right)\mathrm{d}{\mathit{\eta }}_b\right]\\ & =\underset{\mathcal{H}}{\int }\underset{\mathcal{X}}{\int }\left[\left(\underset{\mathrm{0}}{\overset{a}{\int }}\widehat{x}\mathit{\delta }\left({\mathit{\eta }}_a-h^{\prime }\right)\mathrm{d}{\mathit{\eta }}_a\right)\left(\underset{\mathrm{0}}{\overset{b}{\int }}\widehat{x}\mathit{\delta }\left({\mathit{\eta }}_b-h^{\prime }\right)\mathrm{d}{\mathit{\eta }}_b\right)\right.\\ & \left.f_{\widehat{X}|H}\left(\widehat{x}|h^{\prime }\right)f_H\left(h^{\prime }\right)\right]\mathrm{d}\widehat{x}\mathrm{d}{h}^{\prime }\\ & =\underset{\mathrm{0}}{\overset{a}{\int }}\underset{\mathrm{0}}{\overset{b}{\int }}\underset{\mathcal{H}}{\int }\left[\mathit{\delta }\left(h^{\prime }-{\mathit{\eta }}_b\right)\mathit{\delta }\left(h^{\prime }-{\mathit{\eta }}_a\right)f_H\left(h^{\prime }\right)\right.\\ & \left.\mathbb{E}\left({\widehat{X}}^{\mathrm{2}}|H=h^{\prime }\right)\right]\mathrm{d}{h}^{\prime }\mathrm{d}{\mathit{\eta }}_b\mathrm{d}{\mathit{\eta }}_a,\end{array}\end{array}$$

because the domains (0,a), (0,b), and ℋ have the relationship

$$\begin{array}{}\text{(B9)}& (\mathrm{0},a)\subseteq (\mathrm{0},b)\subseteq \mathcal{H},\end{array}$$

it follows that η_a∈ℋ, η_b∈ℋ and $\mathrm{0}<{\mathit{\eta }}_a<{\mathit{\eta }}_b$, so that

$$\begin{array}{}\text{(B10)}& \begin{array}{rl}\mathbb{E}& \left[\widehat{X}\left(a\right)\widehat{X}\left(b\right)\right]=\underset{\mathrm{0}}{\overset{a}{\int }}\underset{\mathrm{0}}{\overset{b}{\int }}\left[\mathit{\delta }\left({\mathit{\eta }}_b-{\mathit{\eta }}_a\right)f_H\left({\mathit{\eta }}_b\right)\right.\\ & \left.\mathbb{E}\left({\widehat{X}}^{\mathrm{2}}|H={\mathit{\eta }}_b\right)\right]\mathrm{d}{\mathit{\eta }}_b\mathrm{d}{\mathit{\eta }}_a\\ & =\underset{\mathrm{0}}{\overset{a}{\int }}{f}_H\left({\mathit{\eta }}_a\right)\mathbb{E}\left({\widehat{X}}^{\mathrm{2}}|H={\mathit{\eta }}_a\right)\mathrm{d}{\mathit{\eta }}_a=\mathbb{E}\left[{\widehat{X}}^{\mathrm{2}}\left(a\right)\right].\end{array}\end{array}$$

Using the algebraic formula for the variance and covariance, we attain by substitution

$$\begin{array}{}\text{(B11)}& \begin{array}{rl}\mathrm{Cov}& \left[\widehat{X}\left(a\right),\widehat{X}\left(b\right)\right]=\mathrm{Var}\left[\widehat{X}\left(a\right)\right]\\ & -\mathbb{E}\left[\widehat{X}\left(a\right)\right]\left(\mathbb{E}\left[\widehat{X}\left(b\right)\right]-\mathbb{E}\left[\widehat{X}\left(a\right)\right]\right),\end{array}\end{array}$$

which is always less than $\mathrm{Var}\left[\widehat{X}\right(a\left)\right]$ unless a=b since the partial VCD is cumulative, yielding larger values at higher altitudes.

Remark

Of particular importance is that these formulae may be used to calculate the expectation and variance values of the partial VCD between two altitudes. The expected value is

$$\begin{array}{}\text{(B12)}& \mathbb{E}\left[\widehat{X}\left(b\right)-\widehat{X}\left(a\right)\right]=\mathbb{E}\left[\widehat{X}\left(b\right)\right]-\mathbb{E}\left[\widehat{X}\left(a\right)\right],\end{array}$$

and the variance is

$$\begin{array}{}\text{(B13)}& \begin{array}{rl}\mathrm{Var}& \left[\widehat{X}\left(b\right)-\widehat{X}\left(a\right)\right]=\mathrm{Var}\left[\widehat{X}\left(b\right)\right]+\mathrm{Var}\left[\widehat{X}\left(a\right)\right]\\ & -\mathrm{Cov}\left[\widehat{X}\left(a\right),\widehat{X}\left(b\right)\right]\\ & =\mathrm{Var}\left[\widehat{X}\left(b\right)\right]+\mathbb{E}\left[\widehat{X}\left(a\right)\right](\mathbb{E}\left[\widehat{X}\left(b\right)\right]\\ & -\mathbb{E}\left[\widehat{X}\left(a\right)\right]).\end{array}\end{array}$$

Appendix C: Bilinear interpolation of background spectrum

To smooth the changes between retrievals in adjacent background cells, we use a bilinear interpolation of the background spectra. For a general quantity (Q), the bilinear interpolation is represented as

$$\begin{array}{}\text{(C1)}& \begin{array}{rl}Q(x,y)& =c_x{c}_{y}Q\left(x_{\mathrm{0}},y_{\mathrm{0}}\right)+\left(\mathrm{1}-c_x\right)c_{y}Q\left(x_{\mathrm{1}},y_{\mathrm{0}}\right)\\ & +c_x\left(\mathrm{1}-c_y\right)Q\left(x_{\mathrm{0}},y_{\mathrm{1}}\right)\\ & +\left(\mathrm{1}-c_x\right)\left(\mathrm{1}-c_y\right)Q\left(x_{\mathrm{1}},y_{\mathrm{1}}\right),\end{array}\end{array}$$

between the corner points (x₀,y₀), (x₁,y₀), (x₀,y₁), and (x₁,y₁), using the scalings $c_x=(x_{\mathrm{1}}-x)/(x_{\mathrm{1}}-x_{\mathrm{0}})$ and $c_y=(y_{\mathrm{1}}-y)/(y_{\mathrm{1}}-y_{\mathrm{0}})$. We interpolate the inverse error covariance matrix S⁻¹ by this formula. However, because our background spectrum is characterized probabilistically as a set of N=10 000 samples (${\mathit{y}}_{\mathrm{bg}}^s\in {\mathrm{\Omega }}_{{\mathit{Y}}_{\mathrm{bg}}}$) of a correlated random vector Y_bg in each seasonal $\mathrm{5}{}^{\circ }\times \mathrm{5}{}^{\circ }$ background grid cell, the samples cannot be interpolated directly by the above formula. Instead we treat ${\mathit{Y}}_{\mathrm{bg}}={\mathit{Y}}_{\mathrm{bg}}(X,Y)$ as a function of a random position where (X,Y) is a discrete random position, taking only the cell corner points as possible values. In this case, X and Y represent longitude and latitude. In particular, we characterize (X,Y) by the probability mass function $p(x_i,y_j)=(-\mathrm{1}{)}^{i+j}(i-c_x\left)\right(j-c_y)$ for $i,j\in \mathit{\{}\mathrm{0},\mathrm{1}\mathit{\}}$, which is simply the corner point weights in the bilinear interpolation formula. Consequently, we generate ${\mathit{Y}}_{\mathrm{bg}}(x,y)={\mathbb{E}}_{(X,Y)}\left[{\mathit{Y}}_{\mathrm{bg}}\right(X,Y\left)\right]={\sum }_{i,j}{\mathit{Y}}_{\mathrm{bg}}(x_i,y_j)p(x_i,y_j)$ by sampling the discrete distribution p(x_i,y_j) to generate the number of samples taken from each of the corner points $n(x_i,y_j)=\left[Np\right(x_i,y_j\left)\right]$, where the bracket represents rounding to the nearest integer. Using this sampling, we generate the samples of Y_bg(x,y) as the collection of each of the n(x_i,y_j) samples from the corners Y_bg(x_i,y_j). This generates a total of N samples for the interpolated background spectrum.

Appendix D: CrIS channels used in SO₂ retrieval

The following CrIS channels are used in this work and are identified by their wavenumber value (cm⁻¹).

D1 For the regular retrieval

1300.0, 1300.625, 1301.25, 1301.875, 1302.5, 1303.125, 1303.75, 1304.375, 1305., 1305.625, 1306.25, 1306.875, 1307.5, 1308.125, 1308.75, 1309.375, 1310., 1310.625, 1311.25, 1311.875, 1312.5, 1313.125, 1313.75, 1314.375, 1315., 1315.625, 1316.25, 1316.875, 1317.5, 1318.125, 1318.75, 1319.375, 1320., 1320.625, 1321.25, 1321.875, 1322.5, 1323.125, 1323.75, 1324.375, 1325., 1325.625, 1326.25, 1326.875, 1327.5, 1328.125, 1328.75, 1329.375, 1330., 1330.625, 1331.25, 1331.875, 1332.5, 1333.125, 1333.75, 1334.375, 1335., 1335.625, 1336.25, 1336.875, 1337.5, 1338.125, 1338.75, 1339.375, 1340., 1340.625, 1341.25, 1341.875, 1342.5, 1343.125, 1343.75, 1344.375, 1345., 1345.625, 1346.25, 1346.875, 1347.5, 1348.125, 1348.75, 1349.375, 1350., 1350.625, 1351.25, 1351.875, 1352.5, 1353.125, 1353.75, 1354.375, 1355., 1355.625, 1356.25, 1356.875, 1357.5, 1358.125, 1358.75, 1359.375, 1360., 1360.625, 1361.25, 1361.875, 1362.5, 1363.125, 1363.75, 1364.375, 1365., 1365.625, 1366.25, 1366.875, 1367.5, 1368.125, 1368.75, 1369.375, 1370., 1370.625, 1371.25, 1371.875, 1372.5, 1373.125, 1373.75, 1374.375, 1375., 1375.625, 1376.25, 1376.875, 1377.5, 1378.125, 1378.75, 1379.375, 1380., 1380.625, 1381.25, 1381.875, 1382.5, 1383.125, 1383.75, 1384.375, 1385., 1385.625, 1386.25, 1386.875, 1387.5, 1388.125, 1388.75, 1389.375, 1390., 1390.625, 1391.25, 1391.875, 1392.5, 1393.125, 1393.75, 1394.375, 1395., 1395.625, 1396.25, 1396.875, 1397.5, 1398.125, 1398.75, 1399.375, 1400., 1400.625, 1401.25, 1401.875, 1402.5, 1403.125, 1403.75, 1404.375, 1405., 1405.625, 1406.25, 1406.875, 1407.5, 1408.125, 1408.75, 1409.375, 1410.0

D2 For the specialized, high-loading retrieval

1300., 1300.625, 1301.25, 1301.875, 1302.5, 1303.125, 1303.75, 1304.375, 1305., 1305.625, 1306.25, 1306.875, 1307.5, 1308.125, 1308.75, 1309.375, 1310., 1310.625, 1311.25, 1311.875, 1312.5, 1313.125, 1313.75, 1314.375, 1315., 1315.625, 1316.25, 1316.875, 1317.5, 1318.125, 1318.75, 1319.375, 1320., 1320.625, 1321.25, 1321.875, 1322.5, 1323.125, 1323.75, 1324.375, 1325., 1325.625, 1326.25, 1326.875, 1327.5, 1328.125, 1328.75, 1329.375, 1330., 1330.625, 1331.25, 1331.875, 1332.5, 1362.5, 1363.125, 1363.75, 1387.5, 1388.125, 1388.75, 1389.375, 1390., 1390.625, 1391.25, 1391.875, 1392.5, 1393.125, 1393.75, 1394.375, 1395., 1395.625, 1396.25, 1396.875, 1397.5, 1398.125, 1398.75, 1399.375, 1400., 1400.625, 1401.25, 1401.875, 1402.5, 1403.125, 1403.75, 1404.375, 1405., 1405.625, 1406.25, 1406.875, 1407.5, 1408.125, 1408.75, 1409.375, 1410.

Appendix E: Probabilistic time series

E1 Probabilistic mass

In general, the total cloud mass can be calculated by integrating the total VCD $\widehat{X}\left(h_{\mathrm{\infty }}\right)$ over the SO₂ cloud region. In the present study, we make this calculation after interpolating the CrIS retrievals onto an equal area grid (grid cells of constant area δA). Because the sets of measurements are normally distributed and assumed independent with means $\mathbb{E}\left[{\widehat{X}}_i\right(h_{\mathrm{\infty }}\left)\right]$ and variances $\mathrm{Var}\left[{\widehat{X}}_i\right(h_{\mathrm{\infty }}\left)\right]$, their sum is also normally distributed (DeGroot and Schervish, 2012):

$$\begin{array}{}\text{(E1)}& M\sim \mathcal{N}\left(\mathbb{E}\left(M\right),\mathrm{Var}\left(M\right)\right),\end{array}$$

where the mean is

$$\begin{array}{}\text{(E2)}& \mathbb{E}\left(M\right)=\mathit{\kappa }\mathit{\delta }A\sum _i\mathbb{E}\left[{\widehat{X}}_i\left(h_{\mathrm{\infty }}\right)\right],\end{array}$$

and variance is

$$\begin{array}{}\text{(E3)}& \mathrm{Var}\left(M\right)=(\mathit{\kappa }\mathit{\delta }A{)}^{\mathrm{2}}\sum _i\mathrm{Var}\left[{\widehat{X}}_i\left(h_{\mathrm{\infty }}\right)\right],\end{array}$$

and M has units of kilotons (kt) of SO₂ and the factor $\mathit{\kappa }=\mathrm{2.8617}\times {\mathrm{10}}^{-\mathrm{11}}$ kt m⁻² DU⁻¹ has been included for dimensional consistency with VCD measured in DU. This parameterizes the total cloud mass as a Gaussian PDF for any period of data coverage.

E2 Probabilistic decay rate coefficient and e-folding time

We treat the above time series of PDFs of SO₂ mass as a random process M_t. As a continuous process, the conversion of SO₂ into sulfur aerosols can be modeled kinetically by the differential equation ${\dot{M}}_t=-k_t{M}_t$, where k_t is the instantaneous decay rate coefficient. Below we generate k_t and the e-folding time ${\mathit{\tau }}_t=k_t^{-\mathrm{1}}$ as random processes from M_t.

To make this calculation in practice, a finite difference formula is needed for ${\dot{M}}_t$. We write this as a general 2α order accuracy central finite difference formula for the first derivative:

$$\begin{array}{}\text{(E4)}& {\dot{M}}_t\approx {\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }t}}\sum _{i=-\mathit{\alpha }}^{\mathit{\alpha }}{\mathit{\delta }}_{t+i}{M}_{t+i},\end{array}$$

where δ_t+i is the central difference scheme coefficient. As before, the weighted sum of normal random variables is also normally distributed. Consequently,

$$\begin{array}{}\text{(E5)}& {\dot{M}}_t\sim \mathcal{N}\left(\mathbb{E}\left({\dot{M}}_t\right),\mathrm{Var}\left({\dot{M}}_t\right)\right),\end{array}$$

where the mean is

$$\begin{array}{}\text{(E6)}& \mathbb{E}\left({\dot{M}}_t\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }t}}\sum _{i=-\mathit{\alpha }}^{\mathit{\alpha }}{\mathit{\delta }}_{t+i}\mathbb{E}\left(M_{t+i}\right),\end{array}$$

and variance is

$$\begin{array}{}\text{(E7)}& \begin{array}{rl}\mathrm{Var}\left({\dot{M}}_t\right)& ={\displaystyle \frac{\mathrm{1}}{(\mathrm{\Delta }t{)}^{\mathrm{2}}}}[\sum _{i=-\mathit{\alpha }}^{\mathit{\alpha }}{\mathit{\delta }}_{t+i}^{\mathrm{2}}\mathrm{Var}\left(M_{t+i}\right)\\ & +\mathrm{2}\sum _{i\ne j}{\mathit{\delta }}_{t+i}{\mathit{\delta }}_{t+j}\mathrm{Cov}\left(M_{t+i},M_{t+j}\right)]\\ & ={\displaystyle \frac{\mathrm{1}}{(\mathrm{\Delta }t{)}^{\mathrm{2}}}}\sum _{i=-\mathit{\alpha }}^{\mathit{\alpha }}{\mathit{\delta }}_{t+i}^{\mathrm{2}}\mathrm{Var}\left(M_{t+i}\right),\end{array}\end{array}$$

where each covariance is zero because each measurement is independent. Because each ${\dot{M}}_t$ is normally distributed, this sequence of means and variances fully parameterizes its random process. To calculate k_t, we also must calculate $\mathrm{Cov}(M_t,{\dot{M}}_t)$:

$$\begin{array}{}\text{(E8)}& \begin{array}{rl}\mathrm{Cov}\left(M_t,{\dot{M}}_t\right)& ={\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }t}}\mathbb{E}\left[M_t\sum _{i=-\mathit{\alpha }}^{\mathit{\alpha }}{\mathit{\delta }}_{t+i}{M}_{t+i}\right]\\ & -{\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }t}}\mathbb{E}\left(M_t\right)\sum _{i=-\mathit{\alpha }}^{\mathit{\alpha }}{\mathit{\delta }}_{t+i}\mathbb{E}\left(M_{t+i}\right)\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }t}}\sum _{i=-\mathit{\alpha }}^{\mathit{\alpha }}{\mathit{\delta }}_{t+i}\mathrm{Cov}\left(M_t,M_{t+i}\right)\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }t}}{\mathit{\delta }}_t\mathrm{Var}\left(M_t\right)=\mathrm{0}\end{array}\end{array}$$

where the last two equalities follow from the independence of each M_t and the fact that the central coefficient (δ_t) in any central finite difference for the first derivative is zero. For non-central differences, this is not zero, the last equality does not hold, and $\mathrm{Cov}(M_t,{\dot{M}}_t)=\frac{\mathrm{1}}{\mathrm{\Delta }t}{\mathit{\delta }}_t\mathrm{Var}\left(M_t\right)$.

With random processes for the mass and mass rate of change calculated we can calculate the decay rate coefficient as a function of these two random processes:

$$\begin{array}{}\text{(E9)}& k_t=-{\displaystyle \frac{{\dot{M}}_t}{M_t}},\end{array}$$

which is a ratio of two Gaussian random processes which may be correlated depending on the finite differencing scheme. The calculation of such a ratio of random variables (Fieller, 1932; Hinkley, 1969) describes the uncertainty of the decay coefficient at each time as a PDF $f_{k_t}\left(k\right)$.

Calculating the PDF for the instantaneous e-folding time ${\mathit{\tau }}_t=k_t^{-\mathrm{1}}$ is performed by applying standard rules for functions of random variables (from DeGroot and Schervish, 2012):

$$\begin{array}{}\text{(E10)}& f_{{\mathit{\tau }}_t}\left(\mathit{\tau }\right)={\displaystyle \frac{\mathrm{1}}{{\mathit{\tau }}^{\mathrm{2}}}}{f}_{k_t}\left({\displaystyle \frac{\mathrm{1}}{\mathit{\tau }}}\right).\end{array}$$

Notably, the distributions are not Gaussian either the decay rate coefficient or the e-folding time .