Articles | Volume 15, issue 20
Research article
20 Oct 2022
Research article |  | 20 Oct 2022

Uncertainty-bounded estimates of ash cloud properties using the ORAC algorithm: application to the 2019 Raikoke eruption

Andrew T. Prata, Roy G. Grainger, Isabelle A. Taylor, Adam C. Povey, Simon R. Proud, and Caroline A. Poulsen

Uncertainty-bounded satellite retrievals of volcanic ash cloud properties such as ash cloud-top height, effective radius, optical depth and mass loading are needed for the robust quantitative assessment required to warn aviation of potential hazards. Moreover, there is an imperative to improve quantitative ash cloud estimation due to the planned move towards quantitative ash concentration forecasts by the Volcanic Ash Advisory Centers. Here we apply the Optimal Retrieval of Aerosol and Cloud (ORAC) algorithm to Advanced Himawari Imager (AHI) measurements of the ash clouds produced by the June 2019 Raikoke (Russia) eruption. The ORAC algorithm uses an optimal estimation technique to consolidate a priori information, satellite measurements and associated uncertainties into uncertainty-bounded estimates of the desired state variables. Using ORAC, we demonstrate several improvements in thermal infrared volcanic ash retrievals applied to broadband imagers. These include an improved treatment of measurement noise, accounting for multi-layer cloud scenarios, distinguishing between heights in the troposphere and stratosphere, and the retrieval of a wider range of effective radii sizes than existing techniques by exploiting information from the 10.4 µm channel. Our results indicate that 0.73 ± 0.40 Tg of very fine ash (radius  15 µm) was injected into the atmosphere during the main eruptive period from 21 June 18:00 UTC to 22 June 10:00 UTC. The total mass of very fine ash decreased from 0.73 to 0.10 Tg over  48 h, with an e-folding time of 20 h. We estimate a distal fine ash mass fraction of 0.73 %± 0.62 % based on the total mass of very fine ash retrieved and the ORAC-derived height–time series. Several distinct ash layers were revealed by the ORAC height retrievals. Generally, ash in the troposphere was composed of larger particles than ash present in the stratosphere. We also find that median ash cloud concentrations fall below peak ash concentration safety limits (< 4 mg m−3) 11–16 h after the eruption begins, if typical ash cloud geometric thicknesses are assumed. The ORAC height retrievals for the near-source plume showed good agreement with GOES-17 side-view height data (R=0.84; bias =0.75 km); however, a larger negative bias was found when comparing ORAC height retrievals for distal ash clouds against Cloud-Aerosol Lidar with Orthogonal Polarisation (CALIOP) measurements (R=0.67; bias =2.67 km). The dataset generated here provides uncertainties at the pixel level for all retrieved variables and could potentially be used for dispersion model validation or be implemented in data assimilation schemes. Future work should focus on improving ash detection, improving height estimation in the stratosphere and exploring the added benefit of visible channels for retrieving effective radius and optical depth in opaque regions of nascent ash plumes.

1 Introduction

Volcanic Ash Advisory Centers (VAACs) require satellite observations to detect and track volcanic ash clouds that pose a threat to aviation. In addition to detection and tracking, VAACs use dispersion models to forecast the position of a volcanic ash cloud. Currently, VAACs are required to provide qualitative, deterministic forecasts indicating the future position of potentially hazardous ash clouds, and satellite detection schemes have been developed in support of this operational requirement (Pavolonis et al.2015a, b). However, according to the International Airways Volcano Watch (IAVW) roadmap, by 2025 (2026/2027 for the full requirement to be made; Chris Lucas, personal communication, 24 May 2022), VAACs will be required to issue quantitative forecasts of ash concentration (ICAO2019). The move to quantitative dispersion model forecasts motivates the need for validation to improve the overall quality of the forecast. The thermal infrared (IR) capabilities of the current generation of geostationary satellite sensors are particularly well suited to this purpose as they offer continuous day and night observations at finer spatial (2 km) and temporal (10 min) resolutions than operational dispersion model output grids (typically 10–25 km and 1–3 h). Thermal IR geostationary satellite observations are useful for quantitative validation (Wilkins et al.2016; Prata et al.2021; Folch et al.2022), source term characterisation (Pouget et al.2013; Van Eaton et al.2016; Prata et al.2020), data insertion (Wilkins et al.2015, 2016; Folch et al.2020; Prata et al.2021), data assimilation (Lu et al.2016; Pardini et al.2020; Zidikheri and Lucas2021; Mingari et al.2022) and source term inversion (Stohl et al.2011; Harvey et al.2020, 2022). In addition, it is important that uncertainties in satellite retrievals are accurately characterised because the VAACs and other users of satellite retrievals increasingly require uncertainty information to correctly interpret, aggregate and utilise the data.

Volcanic ash retrievals from satellite data are possible in the thermal IR because of the high SiO2 content of volcanic ash, which has a strong absorption feature around the 9.5 µm wavelength region (Soda1961; Grainger et al.2013; Prata et al.2019). Volcanic ash clouds can be discriminated from water and ice clouds because the absorption of thermal radiation for ash decreases from 10 to 12 µm, while thermal infrared absorption increases from 10 to 12 µm for water and ice, a property known as the reverse absorption effect (Prata1989a, b). Wen and Rose (1994) and Prata and Grant (2001) demonstrated how estimates of the effective radius and optical depth (at 11 µm) could be obtained on the basis of two split-window brightness temperature measurements centred near 11 and 12 µm. Based on this principle, mass loadings (mass per unit area) can be derived at pixel-scale resolution, and when combined with vertical profile information, ash concentrations can be estimated (Prata and Prata2012).

Uncertainties in satellite-based ash cloud retrieval algorithms are dominated by inaccuracies in the physical model whereby ash cloud properties or state variables (such as cloud-top pressure or height, optical depth and effective radius) are converted into satellite-measured radiances. This is achieved using a radiative transfer forward model (FM) which simulates top-of-the-atmosphere (TOA) radiances based on certain assumptions about the atmospheric state. Numerous simplifying assumptions are needed to evaluate the FM, and each of these assumptions introduce FM uncertainty. Examples include parameterising a cloud layer in terms of its microphysical, optical and geometric properties, assuming whether or not it is plane parallel (i.e. neglecting 3D radiative effects) and parameterising the underlying surface characteristics. Additional uncertainties arise due to uncertainty in ancillary information, often derived from reanalysis data derived from numerical weather prediction (NWP) models, such as profiles of temperature, humidity and trace gases. Uncertainties related to the measurements must also be considered, which include channel noise and co-registration and scene-dependent uncertainties introduced by sub-pixel-scale inhomogeneity.

Uncertainty in the total mass of very fine ash (radius  15 µm) derived using the split-window method has been estimated from sensitivity analyses and is generally found to be 40 %–60 % (Wen and Rose1994; Gu et al.2003; Corradini et al.2008; Prata and Prata2012; Prata et al.2017b). Sources of FM uncertainty that have previously been considered include the assumed ash composition (e.g. Wen and Rose1994; Mackie et al.2014; Prata et al.2019; Deguine et al.2020; Piontek et al.2021b), form of the underlying size distribution (e.g. Wen and Rose1994; Western et al.2015), variation in surface and cloud-top temperature (e.g. Corradini et al.2008; Schneider et al.1999), particle shape (Kylling et al.2014) and meteorological cloud interference (Kylling et al.2015). While broad estimates of uncertainty in the total mass of very fine ash are useful, it is desirable to have pixel-scale, uncertainty-bounded estimates of all retrieved quantities.

Optimal estimation (OE; Rodgers1976, 2000) is particularly useful for solving the noisy inverse problem, where imperfect prior knowledge in the state space is quantified in terms of probability density functions (PDFs), and noise in the measurements is quantified in terms of PDFs in the measurement space. Within the OE framework, one can assume large uncertainties in the priors, such that the solution is influenced mainly by the measurements or set small uncertainties if good quality a priori information is available. Optimal estimation also provides a formalism where state and measurement variables can be easily added or removed, and uncertainties in prior information and the measurements are propagated through the FM equations to estimate uncertainties in all retrieved state variables.

There are several examples of OE retrievals that use thermal-only measurements of volcanic ash, each using different state variables and instrument channels but all aiming to retrieve the ash mass loading and cloud-top height. For example, the Pavolonis et al. (2013) algorithm, used by the Washington, Anchorage and Darwin VAACs, includes the effective cloud temperature, effective cloud emissivity (at 11 µm) and the β ratio (the ratio of effective optical depth at 12 and 11 µm) in the state vector. From these state variables, cloud-top height and mass loading are derived. The measurement vector used by Pavolonis et al. (2013) includes the 11 µm brightness temperature and two brightness temperature differences (BTDs; 11–12 µm and 11–13.3 µm). Uncertainties considered by Pavolonis et al. (2013) include measurement (instrument) error, clear-sky radiance errors (land and water are treated separately) and spatial heterogeneity errors. For the measurement errors in each channel, a fixed, noise-equivalent delta temperature (NEΔT), as reported by the satellite provider, is used. The Francis et al. (2012) algorithm, used by the London VAAC, includes the effective radius, cloud-top pressure and mass loading in their state vector and uses the 11, 12 and 13.3 µm brightness temperatures in their measurement vector. Similar to Pavolonis et al. (2013), Francis et al. (2012) also use fixed NEΔT values to estimate the measurement uncertainty for each channel included in the measurement vector, while FM uncertainty is estimated using statistics derived from long-term, cloud-free satellite radiances. Kylling et al. (2015) present a simple OE scheme that retrieves two state variables (effective radius and optical depth at 11 µm) based on two measurements (11 and 12 µm brightness temperatures) from which mass loading is computed. Cloud-top height is not provided in their paper but could have been inferred from their estimation of cloud-top temperature using a nearby meteorological temperature profile (from a NWP model or sounding). To characterise uncertainty in their retrieval, Kylling et al. (2015) use the combined (FM and measurement) uncertainties provided in Francis et al. (2012).

A common theme amongst existing retrieval schemes is the use of a fixed measurement error (i.e. NEΔT) per channel. However, this noise estimate is only true for the reference temperature given. It is straightforward to allow NEΔT to vary based on the measured brightness temperature, and we highlight this improvement in the present study (Sect. 3.2.2). Another key result from the work of Pavolonis et al. (2013) and Francis et al. (2012) is that the inclusion of the 13.3 µm channel enables good estimates of ash cloud-top height, as shown by validation with the Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation (CALIPSO) satellite lidar measurements presented in their studies. As we will show, the 13.3 µm channel is key for distinguishing between ash clouds in the troposphere and stratosphere. Additionally, no authors have yet published an ash OE scheme that incorporates the 10.4 µm channel in addition to channels centred near 11, 12 and 13.3 µm. The 10.4 µm channel is a window channel less affected by water vapour than the 11 and 12 µm channels (Lindsey et al.2012) and is positioned closer to the 9.5 µm silica absorption band.

We focus on quantifying ash cloud properties and their associated uncertainties using the Optimal Retrieval of Aerosol and Cloud (ORAC) algorithm applied to Advanced Himawari Imager (AHI) observations using the June 2019 Raikoke (Russia) eruption as a case study. The main aim of this study is to provide uncertainty-bounded estimates of optical depth (at 550 nm), effective radius, cloud-top height and mass loading for the Raikoke ash clouds, while simultaneously describing new advances in thermal-only volcanic ash retrievals applied to satellite imager instruments using an OE framework. Specifically, we present a time series of the Raikoke ash cloud properties and discuss our results in the context of existing studies on Raikoke. We present a new method for quantifying ash cloud heights in the troposphere and stratosphere (a common feature of explosive volcanic eruptions such as Raikoke). We also present the first ash retrievals for the Raikoke case which consider multi-layered cloud (ash over water cloud) and present a technique for selecting the best-fitting FM based on measurement cost at solution. We consider uncertainty variation in measured brightness temperatures, which advances previous methods of assuming constant brightness temperature errors. Finally, we explore the advantages of using the 10.4 µm channel in OE retrievals of volcanic ash and demonstrate that its inclusion enables the retrieval of a wider range of effective radii sizes.

2 Data

2.1 Advanced Himawari Imager

The Advanced Himawari Imager (AHI) aboard the Japanese Meteorological Agency's (JMA) Himawari-8 satellite has been in operation since 7 July 2015. The AHI is in geostationary orbit nominally positioned at 140.7 E and completes a full disk scan every 10 min (Bessho et al.2016). The AHI has 16 spectral bands, with spatial resolutions at a nadir of 0.5 km (band 3; 0.64 µm), 1 km (bands 1, 2 and 4; 0.47, 0.51 and 0.86 µm) and 2 km (bands 5–16; 1.6, 2.3, 3.9, 6.2, 6.9, 7.3, 8.6, 9.6, 10.4, 11.2, 12.4 and 13.3 µm). The ORAC algorithm has been applied to level 1b Himawari Standard Data (HSD) files, which contain radiances that have been sampled onto the World Geodetic System 1984 (WGS 84) ellipsoid. The radiances were converted to brightness temperatures using the JMA standard calibration (calibration information block 5 of the header). A subset of the AHI full disk, defined by 135 E–15 W and 40–65 N, has been analysed from 18:00 UTC on 21 June 2019 to 18:00 UTC on 28 June 2019. This corresponds to 1 week of observations from eruption onset.

2.2 Advanced Baseline Imager

To validate ORAC–AHI retrievals of cloud-top height in the near-source volcanic plume, we use the recently published dataset from Horváth et al. (2021a, b), who provide geometric estimations of plume-top height. Specifically, we use the GOES-17 Advanced Baseline Imager (ABI; positioned at 137.2 W) side-view height estimations provided in the Supplement of Horváth et al. (2021b). Uncertainty in these height retrievals is ±500m.

2.3 Cloud-Aerosol Lidar with Orthogonal Polarisation (CALIOP)

To validate the ORAC–AHI retrievals of cloud-top height for distal ash clouds, we use the level 2 lidar products generated from measurements made by the Cloud-Aerosol Lidar with Orthogonal Polarisation (CALIOP) aboard the CALIPSO platform (Winker et al.2009). The level 2, version 4.20, 5 km combined cloud and aerosol layer product (L2_05kmMLay-Standard-V4-20) is used to extract ash cloud heights and geometric thicknesses for validation purposes. The precision with which CALIOP measures the layer top and base height varies with altitude. From −0.5 to 8.2 km the vertical resolution is 30 m and from 8.2 to 20.2 km the vertical resolution is 60 m.

3 Method

3.1 Forward model

The ORAC algorithm is an open-source software initially developed by the University of Oxford and Rutherford Appleton Laboratory (RAL). The Deutscher Wetterdienst (DWD; German Weather Service) has developed the code alongside Oxford and RAL since 2010. The Australian Bureau of Meteorology is also a contributor. Within ORAC, two distinct FM implementations are available to retrieve aerosol or cloud properties (or combined for a joint retrieval). A cloud is considered to be a geometrically infinitesimal layer within the atmosphere, whereas an aerosol layer is considered to be a continuum. Here we use the FM representing cloud because volcanic ash clouds are often observed as well-bounded features in the vertical (e.g. Prata et al.2017a) rather than a well-mixed continuum distributed over a vertical region of the atmosphere. This choice has a practical advantage in that the code for the cloud FM is set up for day and night retrievals, whereas the aerosol FM currently only permits daytime retrievals. Details of the aerosol FM are provided in Thomas et al. (2009), and the cloud FM details can be found in Poulsen et al. (2012) and McGarragh et al. (2018). As we use the ORAC FM for cloud, no further details of the aerosol FM are given here. In the cloud model, the solar and thermal components of terrestrial radiation are considered separately but are minimised simultaneously. As the Raikoke ash clouds dispersed into the atmosphere over several days and nights, we only use thermal channels in our measurement vector to ensure consistency over day and night in the retrievals. As we only use thermal channels in our measurement vector, we focus on the thermal FM within ORAC.

The top-of-the-atmosphere radiance (LTOA), measured by a downward-looking satellite, for a plane-parallel cloud at thermal wavelengths (3–15 µm), can be written as follows (McGarragh et al.2018):

(1) L TOA = L ac + [ L ac R db ( θ v ) + B ( T c ) ϵ ( θ v ) + L bc t db ( θ v ) ] t ac ( θ v ) ,

where Lac is the above-cloud upwelling radiance, Lac is the above-cloud downwelling radiance, Rdb(θv) is the portion of above-cloud downwelling radiance reflected towards the viewing direction of the satellite, θv, B(Tc) is the Planck radiance at the cloud-top temperature, Tc, ϵ(θv) is the cloud emissivity, Lbc is the below-cloud upwelling radiance, tdb(θv) is the upward diffuse transmittance of the cloud, and tac(θv) is the above-cloud transmittance of the atmosphere. Note that, if we assume that the reflected portion of the above-cloud downwelling thermal radiance is negligible (i.e.  LacRdb0) and recognise that tdb(θv)=1-ϵ(θv), then Eq. (1) reduces to the FM formulations used in Pavolonis et al. (2013) and Francis et al. (2012), as follows:

(2) L TOA = ϵ L cld + ( 1 - ϵ ) L clr ,


(3) L clr = L bc t ac + L ac ,


(4) L cld = L ac + B ( T c ) t ac .

Furthermore, if it is assumed that the atmosphere is perfectly transparent and the surface has an emissivity of 1, then Lclr=B(Ts) and Lcld=B(Tc), and we arrive at the original formulation proposed by Prata (1989a) (i.e. LTOAϵB(Tc)+(1-ϵ)B(Ts)). Therefore, although there are significant differences in the practical implementation of the ORAC thermal FM and those of Pavolonis et al. (2013) and Francis et al. (2012), the main difference in its theoretical formulation is the inclusion of non-zero above cloud reflectance of downwelling radiance (as written in Eq. 1). Another important difference is that we include surface temperature, Ts, in the state vector. This means that Lbc must be updated during the retrieval process. For computational efficiency, Lbc is written as a linear expansion in temperature as follows (McGarragh et al.2018):

(5) L bc = L bc,a + ( T s - T s,a ) B ( T s,a ) T s,a ϵ s t bc ( θ v ) ,

where Ts,a is the a priori surface temperature (taken from NWP data), ϵs is the emissivity of the surface, and tbc(θv) is the transmittance from the surface to the cloud layer. The clear-sky radiance and transmittance terms are computed using version 13 of RTTOV (Radiative Transfer for TOVS; Saunders et al.2018). As a pre-processing task, RTTOV is run on atmospheric profiles of temperature, specific humidity and ozone taken from ERA5 reanalysis data (Hersbach et al.2020). The ERA5 data are interpolated in time to match the satellite observation time from a 0.5× 0.5 global grid at 6 h temporal resolution. Clear-sky, above-cloud and below-cloud transmittance and radiance profiles are then passed to ORAC. A cloud layer is inserted into the FM by interpolating the above- and below-cloud radiance and transmittance profiles to the first-guess pressure of the state vector. Since the cloud layer is assumed to be a geometrically (but not optically) infinitesimal layer, this implementation is fast and flexible because, once the pre-processing task is done, any cloud optical properties can be introduced or modified. In addition, this approach allows for both single- and multi-layer FM configurations and delegates the generation of computationally expensive single-scattering cloud properties to offline calculations.

The single-scattering properties for the cloud layer are generated as look-up tables (LUTs) using version 2.1 of DISORT (Stamnes et al.2000). Volcanic ash LUTs are generated as a function of the 550 nm optical depth (τ; 17 grid points from 0–256 with log 10 spacing), effective radius (re; 13 grid points at 0.1 µm, 0.5 µm and 1–15 µm in 1 µm intervals), satellite zenith angle (θv; 10 grid points from 0–90 in 10 intervals) and wavelength (λ; convolved to the relevant channel spectral response function). Ash composition information is accounted for using complex refractive index data taken from the Oxford Aerosol Refractive Index Archive (ARIA,, last access: 24 May 2022). As we do not have refractive index data for the Raikoke ash, we ran the retrieval with three different types of ash (with varying bulk silica content), as reported by Reed et al. (2018), that were sampled from the 2010 Eyjafjallajökull (Iceland) eruption, the 1981 Mount Spurr (Alaska, USA) eruption and the 2008 Chaitén (Chile) eruption. After running the ORAC retrieval for all three ash compositions, we found that the Eyjafjallajökull ash consistently outperformed the other two ash compositions (i.e. lower cost and more retrievals converging). The bulk silica contents for Eyjafjallajökull, Mount Spurr and Chaitén are 58.85, 55.99 and 74.90 wt % (see Prata et al.2019, Table 2). Smirnov et al. (2021) provide bulk silica contents for samples representing the 21–26 June 2019 eruption and show that, for glass compositions of shards from air fall ash, bulk silica contents are mostly between 57–63 wt % (see the total alkali silica diagram in their Fig. 5 top panel). Therefore, the bulk composition for the Eyjafjallajökull ash (58.85 wt %) appears to be consistent with the glass shards of air fall ash reported by Smirnov et al. (2021). We therefore present retrieval results only for the Eyjafjallajökull ash composition. Particles were assumed to be spherical, and the underlying size distribution was assumed to follow a lognormal distribution with a spread of 2. While ash particles are known to be of an irregular shape, here we assume that they are spherical for the following two main reasons: (1) at thermal infrared wavelengths larger than 10 µm, uncertainty due to unknown particle habit (non-sphericity) is expected to have little impact on the retrievals, as discussed by numerous previous authors (Wen and Rose1994; Corradini et al.2008; Clarisse et al.2010; Newman et al.2012; Pavolonis et al.2013; Prata et al.2017b). Yang et al. (2007) show that the impact of non-sphericity at thermal infrared wavelengths is negligible for desert dust (which is similar in many ways to volcanic ash). (2) The exact non-spherical shapes of the Raikoke ash particles under investigation here are unknown. Therefore, approximating their shape by some other irregular shape may introduce further errors than simply assuming a sphere. We recognise that it is possible to find differences between spherical and non-spherical particles if irregular porous objects are compared with spheres. Kylling et al. (2014) found that differences in the total mass uncertainty would increase from 40 % to 50 %. However, it is questionable as to how representative the particle shapes used in the Kylling et al. (2014) study are for the Raikoke ash, and therefore, we cannot conclude that the 10 % uncertainty found by Kylling et al. (2014) would apply here.

ORAC currently allows for the inclusion of two cloud layers in the FM. As there was prevalent stratus cloud during the Raikoke eruption, we ran ORAC in both the single-layer and multi-layer mode. In the multi-layer mode, we tried two FM configurations, namely one with a tightly constrained, low-level (800 hPa) water layer underlying an ash layer and the second with a tightly constrained, mid-level (500 hPa) water cloud underlying an ash layer. We also varied a priori settings for the single- and multi-layer runs (described in Sect. 3.2.3).

3.2 Optimal estimation

The OE technique implemented in ORAC utilises Bayes' theorem so that uncertainties in a priori information can be considered in addition to uncertainties (noise) in the satellite measurements (Rodgers2000). In practice, the goal of OE is to minimise a cost function that is described by a χ2 distribution as follows:

(6) χ 2 = [ y - F ( x , b ) ] T S ϵ - 1 [ y - F ( x , b ) ] + [ x - x a ] T S a - 1 [ x - x a ] ,

where y, x and xa are the measurement, state and a priori state vectors, respectively. F(x,b) is the FM vector (i.e. Eq. 1 converted to brightness temperatures for each satellite channel), which is a function of ancillary information, b, and x. Forward model and measurement uncertainties are contained in the measurement error covariance matrix, Sϵ, and a priori uncertainties are contained in the a priori error covariance matrix, Sa. Here it is worth noting that Sϵ and Sa are assumed diagonal, and so all off-diagonal elements (i.e. the covariances) are zero. Therefore, when making this assumption, one should be careful to select state variables that are independent of each other. To minimise Eq. (6), the ORAC algorithm uses the well-known Levenberg–Marquardt minimisation scheme (Levenberg1944; Marquardt1963). Details regarding the implementation of Levenberg–Marquardt for the present study are provided in McGarragh et al. (2018).

3.2.1 State and measurement vectors

When constructing state and measurement vectors, it is important to consider whether or not the measurements contain enough information about the state that is being retrieved. It is therefore good practice to use at least as many independent measurements as there are state variables. The ORAC state vector contains the following five state variables: optical depth at 550 nm (τ), effective radius (re), cloud-top pressure (pc), surface temperature (Ts) and cloud fraction (f). Experience with using ORAC to retrieve f for meteorological clouds has shown that there is a compensating effect if the optical depth and cloud fraction are simultaneously retrieved. Essentially, the optical depth increases as the cloud fraction is reduced (and vice versa). One approach to addressing this issue is to tightly constrain cloud fraction if it can be estimated from higher-resolution data. Watts et al. (1998) and Poulsen et al. (2012) discuss this in more detail. Investigating this issue would be beyond the scope of the present study, and so for now we assume a cloud fraction of 1 and account for uncertainty due to cloud inhomogeneity (i.e. broken cloud or pixels at cloud edges), in addition to errors due to the plane-parallel cloud assumption, as a forward model error (see Sect. 3.2.2) which stems from the work of Watts et al. (1998). We therefore do not attempt to retrieve f and assume it is always equal to 1. Thus, the state vector used for the ash retrievals presented here contains the following four state variables:

(7) x = log 10 ( τ ) r e p c T s .

Early studies on volcanic ash clouds (e.g. Prata1989a; Wen and Rose1994; Prata and Grant2001) have shown that, for a given Ts and Tc, re varies with the 11–12 µm BTD, and τ varies with the brightness temperature measured at 11 µm (or 12 µm). For opaque ash clouds, the BTD will be ∼0K, and so there will be no information on particle size in this case. Two-channel thermal infrared retrievals work best when the ash cloud is semi-transparent and there is a strong thermal contrast between the surface and the cloud (Prata and Prata2012). The dependence of re on the BTD has also been shown for semi-transparent ice and water clouds (Inoue1985; Yamanouchi et al.1987; Prabhakara et al.1988; Parol et al.1991; Key1995; Cooper et al.2006; Wang et al.2011). This explains why OE schemes attempting to retrieve re and τ from thermal-only measurements often include window channels centred near 11 and 12 µm in the measurement vector. Furthermore, as discussed in Pavolonis et al. (2013) and Francis et al. (2012), the addition of the 13.3 µm to the measurement vector improves the cloud-top pressure (height) estimation, particularly for optically thin upper troposphere/lower stratosphere (UTLS) ash clouds, and the same applies to cirrus clouds (Heidinger et al.2010). Channels centred near 10.4 µm aboard geostationary imagers (e.g. ABI and AHI) have only recently become available, and it is thought that, with the inclusion of the 10.4 µm channel to the measurement vector, additional microphysical information on volcanic ash particles may be extracted (Pavolonis and Sieglaff2012; Pavolonis et al.2020). We therefore include AHI channels 10.4, 11.2, 12.4 and 13.3 µm in our measurement vector so that

(8) y = T 10 T 11 T 12 T 13 ,

where T10, T11, T12 and T13 represent the brightness temperatures in the 10.4, 11.2, 12.4 and 13.3 µm AHI channels, respectively. These channels are advantageous not just for the reasons mentioned above but also because they are unaffected by sulfur dioxide (SO2) absorption, which is not currently included in our FM and was present in abundance for the Raikoke case (Hyman and Pavolonis2020; Prata et al.2021; de Leeuw et al.2021).

3.2.2 Forward model and measurement uncertainty

Forward model uncertainties arise due to assumptions and approximations used to evaluate the FM. Based on previous studies with earlier versions of the ORAC algorithm (Watts et al.1998, 2011; Poulsen et al.2012), we assume fixed uncertainty in the thermal channels of 0.50 K which accounts for uncertainties due to the plane parallel assumption (i.e. 3D radiative effects) and sub-pixel-scale inhomogeneity (Iwabuchi and Hayasaka2002). We also account for errors due to misalignment (co-registration errors) between channels, assuming a fixed uncertainty of 0.15 K for each channel.

To estimate measurement uncertainties in thermal infrared channels, we use the noise-equivalent delta temperature (NEΔT). As mentioned earlier, the NEΔT for satellite imager channels is reported at a particular reference temperature (T0), and thus, it is only accurate for that reference temperature. However, it is straightforward to allow error in the measurements to vary with the measured brightness temperature, which recognises that increased signal (i.e. higher brightness temperatures) will result in reduced noise and vice versa. We therefore compute measurement noise (δTm) for each channel based on the measured brightness temperature per pixel as follows:

(9) δ T m = δ T 0 B ( T 0 ) T 0 / B ( T m ) T m ,

where δT0 is the NEΔT reported by the satellite instrument provider, and B(T0) and B(Tm) are Planck functions evaluated at T0 and Tm. Figure 1 shows how Eq. (9) allows the measurement noise to vary for a range of brightness temperatures for each channel used in our measurement vector. For reference, the NEΔT for each channel is plotted as a horizontal dashed line to demonstrate how uncertainties using a fixed value lead to underestimations of uncertainty for colder brightness temperatures.

Figure 1Relationship between fixed noise uncertainty values (dashed lines) and varying uncertainty noise values (solid lines) derived from Eq. (9). Solid circles indicate reference temperatures (T0) given for the NEΔT estimate.


3.2.3 Lowest cost, a priori settings and first guess

One of the advantages of the ORAC algorithm is that it allows the user to easily modify cloud layer properties for both single-layer and multi-layer FM configurations. Given that there were numerous meteorological clouds (at low- to mid-tropospheric levels) underlying the Raikoke ash clouds, we considered both single-layer ash and multi-layer (ash over water cloud) scenarios. A further consideration was how to deal with local vs. global minima in the cost surface. Considering that the Raikoke ash dispersed into the troposphere and stratosphere (Muser et al.2020; Horváth et al.2021b), and that we use thermal-only channels in our retrieval, the height retrievals are strongly dependent on the temperature profile. Therefore, retrieving heights in the troposphere and stratosphere poses the potential problem of multiple solutions (or multiple minima in the cost surface) due to the inversion of temperature at the tropopause. An additional complication is in the case of an isothermal region in the atmosphere (or a flat cost surface), which is fairly typical of the lower-stratosphere at high latitudes. A nearby radiosonde sounding at around the time of the Raikoke eruption illustrates the problem (Fig. 2).

Figure 2Radiosonde sounding from Kamchatskij station (53.08 N, 158.58 E, 84 m) at 12:00 UTC on 21 June 2019. Data accessed from the University of Wyoming sounding database (, last access: 26 February 2022). ERA5 temperature for the grid box corresponding to Kamchatskij station's location and sounding time is overplotted in blue.


To address issues related to the retrieval of cloud-top pressure, we ran ORAC using five different configurations representing different choices of the a priori pressures in the single- and multi-layer FMs. The a priori pressure settings were chosen based on CALIPSO observations at the beginning of the eruption (see Prata et al.2021, Fig. B3(a)) and are summarised in Table 1. For retrieval configurations with a tropospheric a priori ash layer pressure, the first guess was set to the pressure level where the measured T11 brightness temperature was closest to the ERA5 temperature (searching from the surface to the top of the atmosphere). For ash layers with an a priori pressure level in the stratosphere, the first guess was set equal to the a priori. The a priori settings for τ, re and Ts are summarised in Table 2. The ash layer a priori for τ (at 550 nm) was set to 0.5, which is a typical value for ash cloud retrievals reported in the literature (e.g. Corradini et al.2016). For the ash layer a priori effective radius, we set re to 5 µm, which corresponds roughly to the centre of particle sensitivity for thermal IR channels (Prata and Grant2001). The ash layer a priori uncertainties in τ and re were set to a large number (1×108) to ensure that these parameters were effectively unconstrained by their a priori values. The a priori values for τ and re for water cloud were chosen based on ORAC cloud retrievals (Poulsen et al.2012; McGarragh et al.2018) applied to a stratus deck close to the ash cloud at the beginning of the Raikoke eruption. The water layer a priori uncertainties for pc, τ and re were tightly constrained, such that the measurements only influenced the solution for the upper ash layer in the multi-layer runs.

Table 1A priori and first-guess settings for pc for both multi-layer and single-layer forward model configurations. For the multi-layer configurations, the first number in the a priori column corresponds to the pressure level of ash and the second number to the water cloud. Lower and upper limits within ORAC for pc are 10 and 1200 hPa, respectively.

Download Print Version | Download XLSX

Table 2A priori settings for τ, re and Ts for the ash and water layers considered. All first guesses for these parameters were set to be equal to their a priori values.

Download Print Version | Download XLSX

After all five retrieval configurations were run, the retrieval configuration which resulted in the lowest cost was selected on a per pixel basis to generate the final retrieval product. Figures 3a–e show the cost at the measurement solution at 23:00 UTC on 22 June 2019 for each retrieval configuration on a per pixel basis, and the resulting cost map when the minimum cost from each configuration is selected (Fig. 3f). For this scene, the differences in the five cost values vary over several orders of magnitude in some parts of the ash cloud but can be quite similar (same order of magnitude) in other parts. The most notable differences in cost amongst the five retrieval set-ups are seen when comparing the stratospheric a priori retrieval configurations to the tropospheric a priori retrieval configurations. To illustrate the relative differences in the five cost values, we have generated a forward model flag, where each pixel is coloured according to the retrieval configuration that resulted in the lowest cost. The T11 brightness temperature (Fig. 3g) and the natural colour composite (Fig. 3h) provide contextual information and show that the selected forward model configurations are reasonable. For example, the stratospheric ash over 500 hPa water cloud configuration (light grey pixels) returned the lowest cost for a part of the ash cloud overlying the cold (high) cloud associated with the cyclone.

Figure 3(a)(e) Measurement cost at solution for each of the forward model configurations (annotated above each subplot) used in the present study. (f) Minimum cost per pixel out of the five configurations shown in panels (a)(e). (g) 11.2 µm brightness temperature. (h) Forward model flag (i.e. forward model configuration that resulted in the lowest cost per pixel). The natural colour composite is plotted beneath for context.

In summary, this approach accounts for multi-layer cloud scenarios, multiple local minima in the cost function (in the troposphere and stratosphere) and reduces the impact of flat cost surfaces (isothermal regions) with the use of the a priori uncertainty settings.

3.3 Ash mass loading and uncertainty

To compute ash mass loading (mass per unit area), we use the following standard formulation used by many previous authors (Wen and Rose1994; Prata and Grant2001; Corradini et al.2008; Pavolonis et al.2013):

(10) m l = 4 3 × τ r e ρ Q ext ,

where τ is the optical depth at 550 nm, Qext is the extinction efficiency factor for ash at 550 nm (Qext≈2 at 550 nm), and ρ is the ash particle density, assumed to be 2300 kg m−1. We assumed an ash particle density of 2300 kg m−1 to be consistent with what is used in the Numerical Atmospheric-dispersion Modelling Environment (NAME) model, as an earlier version of the retrievals presented here were used for comparison to the model in Harvey et al. (2022). This ash density was determined for operational use and is therefore a representative average value (Witham et al.2019). We have accounted for uncertainty in this value by allowing for an absolute uncertainty of 300 kg m−1. It is interesting to note that the definition of ash mass loading (Eq. 10) is analogous to the cloud liquid water path found in numerous cloud retrieval studies (see Eq. 24 of Poulsen et al.2012, for example). To compute uncertainty in the ash mass loading, Δml, we assume all variables are independent and sum the error terms in quadrature as follows:

(11) Δ m l m l 2 = Δ τ τ 2 + Δ r e r e 2 + Δ ρ ρ 2 + Δ Q ext Q ext 2 .

Note that all error terms here are retrieved, except for the ash particle density (Δρ) and ΔQext, which is assumed to be negligible compared to the other error terms.

3.4 Ash detection flag

By default, ORAC is run on every level 1b satellite pixel in the full disk image. In some circumstances, retrievals can converge, albeit with a poor fit to the measurements (high cost), even when the FM is not representative of the observation. To avoid these situations and speed up processing times, we only considered pixels within a spatial region from 135 E–15 W and 40–65 N that had a 11–12 µm BTD of less than 0.5 K, which is a fairly loose constraint for ash detection. Pixels were then flagged as “ash” if the water-vapour-corrected BTD, ΔTash, was less than −0.20K. The water vapour correction was applied following the approach of Yu et al. (2002). After this initial ash detection threshold was applied, there were generally two cases that resulted in false positives, i.e. (1) surface inversions and (2) inversions above cloud tops. We removed these cases by setting the pixel to “ash free” if

(12) - 1.25 < Δ T ash < - 0.20  and  T 11 > 275 K ,


(13) - 0.40 < Δ T ash < - 0.20  and  T 11 < 240 K .

These thresholds were chosen based on a manual inspection of the data. A further step to improve the ash flag was to apply an opening morphological 3×3 spatial filter designed to remove isolated pixels unrelated to the ash cloud/plume. Finally, all pixels with satellite zenith angles greater than 75 were ignored, as the plane parallel assumption breaks down at extreme satellite view angles.

3.5 Quality control

To ensure that only the highest quality retrievals were considered for scientific interpretation, we ran a quality control test on each pixel identified by the ash flag (Sect. 3.4). The quality control checks that the retrieval converged all-state variable relative uncertainties were not greater than 100 % and that the retrievals were within a physically sensible range. For the present analysis, the range of valid values considered were 0–15 µm, 0–20 and 0–35 km for effective radius, optical depth and cloud-top height, respectively.

3.6 Gap filling

After quality control, we noticed that gaps in the retrieval fields were appearing where the ash flagging had originally detected ash and that the number of gaps varied with time. In general, the number of gaps increased as the total number of ash-contaminated pixels increased. Given that we have good information (quality-controlled retrievals) adjacent to these gaps, we implemented an algorithm that aims to fill the gaps in the retrieved fields. Specifically, we implemented the Qhull algorithm (Barber et al.1996), which identifies the convex hull of a set of arbitrary points. After finding the convex hull, Delaunay triangulation is used to perform linear barycentric interpolation to fill the missing data. At certain times, the fraction of gap-filled pixels can be significant (reaching as high as  34 %). Figures 4a and b show the ash mass loading at 08:00 UTC on 23 June 2019 before and after gap filling (at this time, the percentage of gap-filled pixels is 23 %). It is important to note that, at the time when the total mass of very fine ash reached its maximum, the fraction of gap-filled pixels was  7 %, meaning that, regardless of gap filling, the maximum total mass estimate is within the uncertainty range estimated here (see Sect. 4.1). The gap-filling algorithm was applied to all retrieved state variables and associated uncertainty fields and the mass loading and uncertainty computed from Eqs. (10) and (11), respectively.

3.7 Parallax correction

Due to the high satellite zenith angles ( 55–75) and high altitude of the volcanic clouds (> 10 km), we also needed to correct the retrievals for parallax. We followed the method of Vicente et al. (2002) to compute the latitude/longitude parallax shift based on the ORAC-retrieved height. We applied the parallax correction to all of the gap-filled, quality-controlled retrievals. The parallax correction resulted in pixel shifts as high as  20 km for observations at the beginning of the eruption (e.g. Fig. 5a). For cases where the parallax shift resulted in two solutions in the same pixel, the set of retrieved fields which corresponded to the higher height was selected. In some instances, the parallax shift can leave behind gaps within the boundaries of the volcanic cloud, meaning that these pixels are being obscured by other parts of the plume. Note that the magnitude of the shift is dependent upon the retrieved height (the shift increases with increasing height for the same viewing angle). To address this, we filled the parallax correction gaps by first applying a 2×2 closing morphological filter to the parallax-corrected ash flag and then filled the retrieval fields using the Qhull algorithm as before (Sect. 3.6).

4 Results and discussion

4.1 Raikoke source term and long-range ash transport

The June 2019 eruption from the Raikoke volcano (48.292 N, 153.25 E, 551 m) is described in detail by McKee et al. (2021). Here we provide an overview of the volcanic ash emissions based on AHI satellite measurements and the ORAC retrieval results. According to AHI measurements, the Raikoke eruption began at around 18:00 UTC on 21 June 2019. What followed was a series of explosive eruptions characterised by sharp decreases in brightness temperatures over the volcano. Nine explosive eruptions can be clearly identified in the AHI time series data, with several smaller events more easily identified with the aid of true colour and thermal imagery.

Figure 4(a) Mass loading retrieval before the gap-filling processing step. (b) Mass loading after the gap-filling processing step.

Figure 5 shows a high-resolution (10 min) time series of the ORAC cloud-top height retrievals for the Raikoke eruption sequence where we have relaxed the BTD thresholds (i.e. removed the threshold conditions in Eqs. 12 and 13) to retrieve the height in the opaque parts of the plume (Fig. 5a). Note that the effective radius retrievals within the opaque regions of the cloud are not reliable when using thermal-only measurements, as BTDs close to zero mean that the solution space cannot be interpolated (i.e. there is no information on particle size for opaque plumes). The time series shows maximum heights (and associated uncertainty) within a search radius of 7.5 km from the volcano (Fig. 5b). The parallax shift is large ( 20 km) for the initial Raikoke plume (grey-shaded regions), demonstrating that, without a parallax correction, significant errors could be introduced when comparing height and location (latitude/longitude) to other satellite datasets or dispersion model output. To compare the ORAC cloud-top height retrievals to the GOES-17 side-view heights, we varied the search radius until the optimal 7.5 km radius was found. We define optimal as the search radius that resulted in the closest match (minimised sum of squared differences) to the GOES-17 side-view height data. This approach ensures a robust comparison with the GOES-17 data, as the exact coordinates of the side-view heights are not provided in the Supplement of Horváth et al. (2021b). The first six pulses were short-lived (eruption duration from 10–40 min), with the seventh being the largest continuous ash emission, lasting almost 4 h. Following a pause of 1 h, two large explosive eruptions occurred at around 03:40 and 05:30 UTC on 22 June 2019.

The maximum height retrieved during the explosive phase of the eruption was 14.2 km (above sea level – a.s.l.), with numerous pulses injecting ash into the stratosphere (tropopause was typically 11.2 km; Fig. 5b). The major phase of the eruption had subsided by 10:00 UTC on 22 June. Based on the ORAC time series, we estimate an eruption duration of 11.7 h and a median plume-top height of 10.7 ± 1.2 km (a.s.l.).

The GOES-17 side-view retrievals, with an estimated uncertainty of ± 500 m (Horváth et al.2021a), serve as validation for heights estimated with ORAC for the initial eruption. The ORAC heights mostly agree well with the GOES-17 side-view heights, with some notable disagreement during the largest duration eruption and the final two pulses (Fig.5b). These differences are not surprising, as the heights reported at these times using the GOES-17 side-view method correspond to small-scale turrets overshooting the main umbrella (see Supplement of Horváth et al.2021b, for details). While ORAC does allow for overshoots by extrapolating the lapse rate (determined from two levels just beneath the tropopause) into the stratosphere (McGarragh et al.2018), the retrieval is relying on thermal infrared measurements that are coarser in spatial resolution compared to the visible channels used by the GOES-17 side-view height analysis. The ORAC heights are, however, an improvement to simply matching a brightness temperature to a NWP profile. Figure 5b shows the cloud-top heights derived using the minimum brightness temperature method described in McKee et al. (2021) (hereafter BTmin method). Note that McKee et al. (2021) used ERA-Interim data, whereas we use ERA5 here. The ORAC and BTmin methods give quite similar results when compared to the GOES-17 heights (similar correlation coefficient and precision). However, there is a significant bias of −1.92km for the BTmin method, while the ORAC cloud-top height bias compared to GOES-17 is −0.75km. This result is to be expected because the BTmin method will, by definition, never produce heights in the stratosphere. While ORAC effectively matches the height to the temperature profile, it accounts for other factors (e.g. cloud emissivity/transmission, viewing angle, above cloud transmission and information from the 13.3 µm channel) compared to matching the raw 11 µm brightness temperature. Furthermore, ambiguity associated with multiple height solutions (due to temperature inversions) is accounted for using optimal estimation by selecting the lowest cost from the output of the five ORAC model run configurations (described in Sect. 3.2.3).

Figure 5(a) Parallax-corrected ORAC ash cloud-top height retrievals at 02:00 UTC on 22 June 2019. Grey shaded regions indicate the parallax shift. The blue circle indicates a 7.5 km radius around the volcano (red triangle). A black dashed box indicates the search region used to find minimum brightness temperatures in McKee et al. (2021). (b) A high temporal resolution (10 min) time series of ORAC ash cloud-top height retrievals for the Raikoke eruption. The ash cloud-top height for the ORAC retrievals (blue dots) in the time series corresponds to the maximum height within at 7.5 km radius of the volcano. Uncertainties associated with the ORAC heights are indicated as light blue shading around the blue dots. Black dots indicate cloud-top heights determined by matching the minimum brightness temperature within the black dashed bounding box in panel (a) to the closest ERA5 grid box profile to the volcano. ERA5 temperature profiles were linearly interpolated vertically and in time to match the minimum AHI brightness temperature every 10 min (i.e. BTmin method). Orange dots indicate the GOES-17 side-view heights taken from Horváth et al. (2021b). Orange error bars indicate ± 500 m. A red dotted line indicates the lapse rate tropopause.

Figure 6ORAC retrievals for (a) mass loading, (b) ash cloud-top height, (c) ash effective radius and (d) ash optical depth at 550 nm. Observation times are annotated at the top of each column in the plot. The blue line on each plot indicates the AHI satellite zenith angle of 75.

Figure 6 shows four scenes (observation times) that illustrate the time evolution and long-range transport of the Raikoke ash cloud. The ORAC mass loading (Fig. 6a) and cloud-top height (Fig. 6b) retrievals show that ash dispersed primarily to the east, with a tropospheric branch of the ash cloud separating toward the south and a stratospheric portion of the ash cloud wrapping up in a cyclone that eventually resulted in ash transport toward the north. The ORAC retrievals show that the tropospheric branch of the ash cloud generally had higher mass loadings (10–30 g m−2) and effective radii (4–7 µm), whereas the stratospheric portion had comparatively lower mass loadings (1–10 g m−2) and effective radii (0.5–2.5 µm). The tropospheric ash cloud maintained heights of 5–8 km, whereas stratospheric ash cloud-top heights remained just above the tropopause at 11–13 km. At 08:00 UTC on 23 June (Fig. 6b second column), there are three distinct levels in the ash cloud, i.e. a mid-tropospheric portion (3.5–4 km), a longitudinally extended, lower-stratospheric region (12 km) and an upper-tropospheric region (6–8 km). By 08:00 UTC on 24 June the mid-tropospheric portion of the ash cloud is no longer detectable with the ash flag, while the upper-tropospheric and lower-stratospheric regions continued to disperse toward the north, reaching extreme satellite view angles (where the ORAC retrievals are not possible) by 08:00 UTC on 25 June.

Figure 7a shows the time series of the total mass of very fine ash from 21 June at 18:00 UTC to 25 June at 18:00 UTC. Although our analysis period covers 7 d, we found that the retrievals began to fail or gave spurious results after 4 d, suggesting that median 550 nm optical depths of  0.1 (Fig. 7b) are at or close to the detection limit for the ORAC ash retrievals applied to AHI. The peak of the total mass time series was reached  13 h after the eruption began (at 07:00 UTC on 22 June 2019), with a total mass of 0.73 ± 0.40 Tg. After reaching its maximum, the total mass decreases with an e-folding time of 20 h. From 24 June onward, approximately 0.1 Tg of the ash remained in the atmosphere, based on the present ash detection and retrieval scheme. The ORAC estimate of the total mass of very fine ash is somewhat lower than existing estimates that also use AHI data to retrieve the total mass (e.g. 1.1 ± 0.7 Tg in Muser et al.2020). One reason for this is due to the difference in ash particle density assumed here (2300 ± 300 kg m−3 and the 2600 kg m−3 value assumed in Muser et al. (2020). Other differences include differing refractive index data, size distribution assumptions, ash detection thresholds, the assumption of a single layer of ash vs. multi-layer cloud/ash scenarios and the use of two channels to retrieve optical depth and effective radius compared to the four channels used in the present study.

The time series of median effective radius is shown in Fig. 7c and reveals that larger particles ( 6 µm) dominated the plume during the first  12 h, with a transition to smaller particles ( 3 µm) from 12–24 h after eruption. After 60 h post-eruption, there is an apparent increase in median effective radius from 3–4 µm. However, given that the median satellite zenith angles (Fig. 7c right axis) were > 70 at this time, it is likely that this increase is due to a retrieval artefact rather than a real increase in particle size in the volcanic cloud. Gu et al. (2005) also found that geostationary satellite retrievals at zenith angles of  70 resulted in an overestimation of effective radius when compared with close-to-nadir (Moderate Resolution Imaging Spectroradiometer, MODIS, and Advanced Very High Resolution Radiometer, AVHRR) retrievals. Factors contributing to retrieval artefacts include the violation of the plane parallel cloud assumption, 3D radiative transfer effects, limb darkening and pixel distortion due to extreme viewing geometry. Parts of the ash cloud may also be exiting the AHI field of view at this time, which could also lead to an abrupt change in the median effective radius.

The time series for median cloud-top heights in the troposphere and stratosphere are shown in Fig. 7d. We chose to plot the median cloud-top heights separately for the troposphere and stratosphere due to the distinct levels of ash that formed following the eruption. Cloud-top heights in the troposphere generally increased from 5–8 km in the first 12 h, followed by a decrease to 5 km after 24 h of atmospheric residence. After this period, much of the ash had fallen out or was undetectable (see Fig. 7a), with the median tropospheric heights varying from 5–7 km (a.s.l.). The median stratospheric cloud-top heights remained fairly constant at 12 km throughout our analysis period, with some minor variation (11–13 km). These heights are close to the chosen a priori stratospheric height (200 hPa) and probably reflect the fact that the measurements had little influence on the retrieved solution due to the isothermal nature of the lower-stratosphere (Fig. 2).

Figure 7Time evolution of (a) total very fine ash mass, (b) median ash optical depth at 550 nm, (c) median ash effective radius and (d) median ash cloud-top height for the first 4 d of atmospheric residence following the Raikoke eruption. The shaded region in panel (a) represents the total mass uncertainty computed from the mass loading uncertainty fields (Eq. 11), and shaded regions in panels (b), (c) and (d) represent the median of the uncertainty in τ, re and hc, respectively. The right axis of panel (c) shows the time series of the median satellite zenith angle (bright blue dashed line), and the shaded region indicates times where the median zenith angle exceeds 70.


4.2 Total mass erupted and distal fine ash mass fraction

The total mass of very fine ash (radius  15 µm) is an important piece of information for dispersion modellers attempting to forecast long-range, fine ash transport. The London VAAC parameterises the source term using an estimate of the plume-top height above vent level (H) converted to a mass eruption rate (M˙), using the empirical fit determined by Mastin et al. (2009). The relationship has seen wide usage in the scientific literature due to its simplicity and the fact that it only requires one input readily determined from satellite data. This power law relationship relates the total mass erupted (all particle sizes) to the maximum plume height (above vent level), and therefore, the distal fine ash mass fraction needs to be set in order to simulate very fine ash transport and dispersion. Typical distal fine ash mass fractions used by the London VAAC range from 1 %–5 % (Dacre et al.2013). Based on a comparison between well-constrained values of M˙, eruption duration and satellite retrievals, Gouhier et al. (2019) found that the distal fine ash mass fraction varies by  2 orders of magnitude (0.1 %–6.9 %) and decreases with increasing M˙. At the wavelengths used here, the retrievals effectively measure the volume of the volcanic ash. For the distal plume, this signal is dominated by fine ash effective radii in the range 0.5–9 µm. For the Raikoke eruption, Osborne et al. (2022) report a total mass of 300 Tg, based on the Mastin et al. (2009) relationship (M˙=140.8H1/0.241), assuming a plume-top height of 15 km and continuous emission from 21:00 UTC on 21 June to 03:00 UTC on 22 June (eruption duration of 6 h). It is not clear how they arrived at the 300 Tg figure, as a constant plume height of 15 km (a.s.l.; 14.45 km above vent level) for a 6 h duration equates to a total mass of 198 Tg when using Mastin et al. (2009)'s empirical fit. Osborne et al. (2022) caution that their total mass estimate should only be taken as a representative figure that is likely an upper bound on the total mass due to the complexity of the eruption source (i.e. numerous pulses, pauses and variations in plume height). The GOES-17 side-view time series data show that the plume heights varied with time from 9–14 km, with overshoots reaching 15–16 km. It is likely that using a constant height of 15 km (a.s.l.) will result in a significant overestimate of the total mass, especially considering the power law relationship between M˙ and H. Despite this, the Osborne et al. (2022) estimate is closer to the lower range of McKee et al. (2021), who determined total masses ranging from 287–672 Tg (average value of 439 Tg) based on more sophisticated plume modelling constrained with plume heights ranging from 10–12 km (a.s.l.). A possible reason for this discrepancy may be due to the fact that Osborne et al. (2022) have underestimated the eruption duration. They assume that the Raikoke eruption started at 21:00 UTC on 21 June and ended at 03:00 UTC on 22 June, which neglects five significant eruptions at around 18:00, 18:50 and 19:40 UTC on 21 June and at around 03:40 and 05:30 UTC on 22 June (see Fig. 5b). McKee et al. (2021) derive a more detailed eruption sequence of events, showing eruptive activity (recorded by infrasound, lightning and AHI data) from around 18:00 UTC 21 June to 10:00 UTC on 22 June. However, those authors note that it is possible that they have systematically underestimated the plume-top height for the Raikoke eruption. They select the coldest pixel from T11 and match it to the temperature corresponding to the highest height below the tropopause (using ERA-Interim data), despite multiple lines of satellite-based evidence showing ash in the stratosphere at the beginning of the eruption. Therefore, the total mass estimates from Osborne et al. (2022) and McKee et al. (2021) broadly agree, but this could be because of compensating errors. Osborne et al. (2022) overestimate the plume height and underestimate the source duration, whereas McKee et al. (2021) underestimate plume-top height.

To understand how the ORAC retrievals translate and compare to existing estimates of the total mass, we converted the time series of ORAC heights (Fig. 8b) to a time series of M˙, based on the Mastin et al. (2009) relationship. We then estimated the total mass by integrating over the M˙ time series. After propagating the ORAC-derived uncertainties in height through the Mastin et al. (2009) equation (see Appendix A), we obtain an estimate of 101 ± 67 Tg for the total mass erupted. Comparing this figure to the maximum total mass of very fine ash derived from the ORAC mass loading retrievals, we estimate a distal fine ash mass fraction of 0.73 %± 0.62 %. Distal fine ash fractions of this magnitude are lower than what is currently used by the London VAAC (1 %–5 %) and adds support to the Gouhier et al. (2019) finding that distal fine ash mass fractions should be set depending on the eruption style.

Figure 8Time evolution of median ash concentrations produced by the Raikoke eruption, assuming different geometric thicknesses. The ICAO peak ash concentration safety limits at 2 and 4 mg m−3 are indicated by dashed horizontal lines. The 0.2 mg m−3 horizonal dashed line indicates the Prata and Prata (2012) detection limit (0.2 g m−2) converted to an ash concentration, assuming a 1 km geometric thickness.


Figure 9(a) Validation results for the ORAC height retrievals compared to the height validation data derived from GOES-17 (side-view heights). (b) Same as panel (a) but for CALIOP. The heights for CALIOP represent mid-layer heights for the top layer detected by the MLay product (with the error bars representing the top and base heights).


4.3 Ash mass loadings vs. ash concentrations

Volcanic ash mass loadings represent the column mass per unit area and are distinct from ash concentrations (mass per unit volume). However, while passive imager measurements do not resolve the ash layers vertically, the mass loadings can be converted to ash concentrations if a geometric thickness is assumed or measured (Sears et al.2013). Prata and Prata (2012) demonstrated that ash concentrations could be derived from ash mass loading retrievals when combined with vertically resolved measurements from CALIOP. While only a few intersections have been identified for CALIOP for the Raikoke eruption, previous studies provide an indication of typical ash cloud geometric thicknesses. Prata et al. (2015) found that ash clouds produced by the 2008 Chaitén eruption were  0.3–0.7 km thick, according to CALIOP observations. Winker et al. (2012) presented CALIOP observations of the 2012 Eyjafjallajökull ash clouds and derived geometric thicknesses of 0.4–1 km. Prata et al. (2017a) studied both ash and sulfates and demonstrated that geometric thicknesses for the ash clouds produced by the 2011 Puyehue-Cordón Caulle eruption were 1.82 ± 0.55 km. Based on comparisons made with CALIPSO in the present study (Sect. 4.4), the Raikoke ash cloud geometric thicknesses were 1.04 ± 0.56 km. Figure 8 shows the median ash concentrations that would be derived if geometric thicknesses of 0.5, 1 and 2 km were assumed. The time series demonstrates that, after 16 h, the median ash concentration (for all geometric thicknesses assumed) would be below what ICAO regards as a high ash concentration level (4 mg m−3).

Figure 10(a) Parallax-corrected ORAC ash cloud-top height retrievals at 01:30 UTC on 23 June 2019. The grey shaded region indicates the parallax shift. The black line indicates the CALIOP ground track, with green indicating the section plotted in panel (c). (b) Correlation between ORAC ash cloud-top heights and CALIOP mid-layer heights. The colours of the data points represent the feature type identified by the CALIOP vertical feature mask (blue for cloud, orange for tropospheric aerosol and brown for stratospheric aerosol). (c) CALIOP level 2 MLay feature classification flags with collocated ORAC heights overplotted.

4.4 Cloud-top height validation

Figure 9a and b show the validation results for the ORAC ash cloud-top heights compared against GOES-17 and CALIOP. In total, we found 115 collocations between the ORAC heights and the validation data derived from CALIOP and GOES-17. There is generally good agreement between the ORAC heights and the GOES-17 data for heights retrieved in the near-source plume at the beginning of the eruption (R=0.84, bias =−0.75km). However, there is a notable negative bias (−2.67km) and poorer agreement (R=0.67) for the ORAC heights retrieved in the distal ash clouds when compared against CALIOP. In both comparisons (Fig. 9a and b), ORAC showed limited skill in retrieving heights in the stratosphere. Negative biases between thermal IR height retrievals and lidar-derived heights have been observed before (Pavolonis et al.2013; Francis et al.2012) and is explained by the fact that the effective thermal emission height of the cloud is generally lower than the cloud-top detectable by lidar backscatter measurements. We tried to account for this by using the CALIOP mid-layer heights (rather than layer-top heights) from the uppermost layer in the MLay product; however, the negative bias persisted. The precision (standard deviation of the difference between the ORAC and validation heights) for the near-source plume heights (GOES-17 comparison) is 1.78 km, which is comparable to, but somewhat higher than, existing thermal IR height retrieval schemes (cf. 1.48–1.64 kmPavolonis et al.2013). However, a key difference here is that we are validating height retrievals in the stratosphere and troposphere, whereas previous validation studies only considered height retrievals in the troposphere. The precision of ORAC heights for the distal ash (CALIOP comparison) was higher (2.75 km) than the precision of ORAC heights retrieved for the near-source plume.

Figure 10 shows a CALIPSO overpass where we identified the largest number of height collocations with the ORAC height retrievals. This observation also serves as an important test case for height retrievals in the troposphere and stratosphere. In general, the height retrievals in the troposphere are underestimated (reflecting the negative bias seen in Fig.9b). The stratospheric height retrievals show very good agreement with the CALIOP observations; however, some ORAC height retrievals returned heights in the mid-troposphere (5–6 km) when there was a stratospheric feature in the CALIOP data from 12–13 km. For these cases, the CALIOP feature optical depth at 532 nm for the top layer was less than 0.05, and so it is likely that the stratospheric volcanic aerosol was too optically thin for the OE to determine the correct height based on the thermal radiances measured by the window channels. Additionally, if the ORAC height retrievals are underestimated, the parallax shift will also be underestimated (for the same satellite zenith angle), leading to collocation errors between AHI and CALIOP. The cloud-top height and associated uncertainty for retrievals in the stratosphere did not deviate significantly from their a priori values. This result means that the measurements had little influence on the stratospheric height retrievals. Although the a priori height in the stratosphere (200 hPa) was chosen based on CALIOP observations at the beginning of the eruption, this result highlights both the difficulty in determining height in an isothermal lower stratosphere and the value in setting a representative a priori in such conditions.

Figure 11(a) Clear-sky transmittance profiles for the T10, T11, T12 and T13 AHI channels taken from the ORAC/RTTOV pre-processor output at 18:00 UTC on 22 June 2019. (b) Same as panel (a) but for the weighting functions. (c) ERA5 temperature profile at the same location as the clear-sky transmittance profiles. (d) Same as panel (c) but for specific humidity.


4.5 Importance of the 13.3 µm channel

The Raikoke case study raises several challenges for reliably retrieving ash cloud-top heights in the troposphere and stratosphere from thermal-only satellite measurements. In particular, there are two main issues, i.e. (1) the double solution due to the inversion of temperature at the tropopause and (2) an isothermal lower stratosphere. We addressed these issues by selecting the lowest cost from a set of FM configurations that used a priori pressures in the troposphere and stratosphere (described in Sect. 3.2.3). To explore the robustness of this method, we ran the five FM configurations again but removed T13 from the measurement vector. We found that, when T13 was not included, the ability to select the correct retrieval solution was lost. Essentially, for retrievals without the T13 channel, the cost for tropospheric and stratospheric height retrievals was very similar, but the stratospheric height solutions always returned the lowest cost.

Given that the information used to retrieve height comes from the thermal channels supplied to the measurement vector, the retrieval is dependent on their respective weighting functions. Figures 11a and b show the transmittance profiles and weighting functions for a clear atmosphere for each of the channels used in the measurement vector. Comparing Fig. 11b with d shows that T12 is most affected by water vapour in the lower troposphere (p>500hPa), followed by T11 and T13, with T10 being the least affected by water vapour. The T10, T11 and T12 channel weighting functions go to zero at  300 hPa, while the T13 band follows the decrease in the temperature profile (Fig. 11c) and remains non-zero up to  10 hPa. The T13 weighting function follows the temperature profile because it is sensitive to CO2 absorption, which is well-mixed in the atmosphere. Crucially, this variation in the weighing function in the T13 channel from 300 to 100 hPa (and above) is what allows the retrieval to distinguish between a cloud layer placed in the troposphere vs. the stratosphere based on cost. Without information from the T13 channel, there is very little difference in TOA radiance for the simulated T10, T11 and T12 channels for a cloud placed at 500 hPa (troposphere) vs. 200 hPa (stratosphere), and so information from the T13 is key for distinguishing between these cases. The weighting functions also show that, because the lower stratosphere is isothermal, the difference in the T13 TOA radiance for a cloud layer placed at 200, 100 and 50 hPa would be very small, meaning that, for a given measurement, the cost would be very similar in each case. This also explains why the retrieved heights in the stratosphere do not deviate significantly from the chosen a priori at 200 hPa (Fig. 10c).

4.6 Importance of the 10.4 µm channel

Figure 12 illustrates the importance of including the 10.4 µm channel in the measurement vector for volcanic ash retrievals using thermal-only channels. By comparing Fig. 12a and b with c, one can see that the T11T12 and T10T11 BTDs are sensitive to different effective radii sizes. The ORAC effective radius retrievals show that small particles are prevalent in the northern part of the volcanic cloud, while larger particles are prevalent in the southern part. Inspection of Fig. 12d shows that, if we were to use a measurement vector that does not include the T10 channel, then the solution space would be restricted to effective particle sizes larger than  2 µm and smaller than  7 µm for this particular scene. In addition, without the T10 channel, effective radii solutions would be found for particles in the northern part of the ash cloud but would be comparatively larger, leading to higher estimates of the total mass. On the other hand, Fig. 12e shows that, if T12 was left out of the measurement vector, and T10 was included, then the retrieval would be sensitive to the smaller particles (for this scene  1–2 µm) but restricted to particle sizes smaller than  5 µm. In addition, a negative T10T11 BTD would not detect the southern part of the ash cloud (Fig. 12b). Therefore, by including T10, T11 and T12 in the measurement vector (Fig. 12f; Supplement), we are able to exploit the particle size information in both the T11T12 and T10T11 BTDs to optimally retrieve a wider range of effective radii sizes than existing retrieval algorithms that only exploit two channels for particle size information (e.g. Corradini et al.2008; Prata and Prata2012; Francis et al.2012; Pavolonis et al.2013).

Figure 12(a) Brightness temperature difference between T11 and T12 at 12:00 UTC on 22 June 2019. The red triangle indicates the location of Raikoke. (b) Same as panel (a) but for a difference between T10 and T11. (c) ORAC effective radius retrieval. (d) Two-dimensional ORAC look-up table for the complex refractive index of the 2010 Eyjafjallajökull ash (Reed et al.2018) plotted in the T11 and T11T12 solution space. The grey data points correspond to pixels within the dashed line bounding box annotated in panel (a). The lines of the constant effective radius are plotted as blue solid lines, and the lines of constant optical depth (at 550 nm) are plotted as dashed orange lines. (e) Same as panel (d) but for the T11 and T10T11 solution space. The grey data points correspond to the bounding box in panel (b). (f) Same as panels (d) and (e) but for the three-dimensional T11, T11T12 and T10T11 solution space.


To understand why the combination of 10.4, 11.2 and 12.4 µm leads to the ability to retrieve a wider range of particle sizes compared to two-channel techniques, it is instructive to consider the heuristic model proposed by Prata and Grant (2001) that explains the relationship between the volume extinction coefficient, effective radius and the BTD for a two-channel retrieval. Essentially, the particle size information is captured by the ratio (β) of volume extinction coefficients at two different wavelengths (k1 and k2) within the thermal infrared window (β=k2/k1). In general, when β>1, the BTD is positive, indicating ice or water. If β=1, then the BTD is 0, and we have no information on particle size. If β<1, then the BTD is negative, and we expect volcanic ash particles. Figure 13 shows how β varies with effective radius for different channel combinations for the LUTs used in the present study (generated assuming a lognormal size distribution, spherical particles and Eyjafjallajökull ash complex refractive index). The Mie calculations show that the extinction coefficient ratio of the 11.2 and 12.4 µm channels contains information on particles sizes from 0.1–9 µm effective radius; however, there is ambiguity (multiple solutions) once the effective radius reaches 1 µm. In other words, the BTD can be the same for effective radii in the range from 0.1 to 1 µm and  1 to 3.5 µm. This is not the case for the 10.4 and 11.2 µm combination, where β continues to decrease as the effective radius reaches 0.5 µm before increasing again. Similar behaviour is seen for the 10.4 and 12.4 µm channel combination. In addition, the channel combination of 10.4 and 11.2 µm will only return negative BTDs (β<1) for effective radii up to 5 µm. This prediction of the Mie theory explains why only the northern part of the Raikoke plume is detected by the 10.4–11.2 µm BTD, whereas the whole plume is detected by the 11.2–12.4 µm BTD (compare Fig. 12a and b). Overall, the combination of 10.4, 11.2 and 12.4 µm allows for an unambiguous retrieval of the effective radius from 0.5–9 µm.

5 Conclusions

In this study, we have presented uncertainty-bounded estimates of volcanic ash cloud-top height, optical depth at 550 nm, effective radius and mass loading for the June 2019 Raikoke eruption. We found the Raikoke eruption injected 0.73 ± 0.40 Tg of very fine ash into the troposphere and stratosphere. After reaching its maximum,  90 % of the total mass was removed from the atmosphere over 48 h (e-folding time of 20 h), with 0.10 Tg detectable in the atmosphere for at least 4 d, corresponding to median ash concentrations of  0.2–2 mg m−3 (depending on the geometric thickness assumed). The distal fine ash mass fraction was estimated to be 0.73 ± 0.62 % based on the total very fine ash mass retrieved and the ORAC cloud-top heights converted to mass eruption rates based on the Mastin et al. (2009) relationship. Our analysis shows that the Raikoke source term is highly complex, meaning that eruption source parameters (i.e. injection height, duration and mass eruption rate) must be carefully considered when attempting to model this eruption. For example, if a continuous eruption source with a constant maximum height is assumed, then it is likely that the total mass will be overestimated due to the numerous pauses and variations in height during this eruption sequence (e.g. Bruckert et al.2022). Even if the duration of the eruptions is accurately captured, underestimates or overestimates could occur if the plume height is not allowed to vary between the troposphere and stratosphere with time.

The ORAC algorithm represents several advances in thermal IR-based ash retrieval algorithms applied to geostationary satellite measurements. Advances include a better characterisation of measurement noise that is allowed to vary with the measured brightness temperature, the ability to distinguish between heights in the troposphere and stratosphere based on cost (with the inclusion of the T13 channel), the retrieval of a wider range of effective radii sizes (with the inclusion of the T10 channel) and accounting for underlying meteorological clouds in the FM. The ash cloud-top height retrievals representing the near-source plume showed good agreement (R=0.84) when compared against GOES-17 side-view height data but showed a notable negative bias (−2.67km) for the distal ash clouds when compared against CALIOP data. Caution must be exercised when interpreting the ORAC ash cloud-top heights in the stratosphere for the distal ash clouds, as the retrieved solutions deviated very little from their a priori values due to the isothermal nature of lower stratosphere. One improvement that could be made would be to use additional information, such as from dispersion model simulations, to set tightly constrained a priori pressure fields where ash is detected by AHI. This approach would essentially be using information from the wind fields in the stratosphere to overcome a flat cost surface due to the isothermal lower stratosphere. Other improvements include using visible channels in the measurement vector to determine optical depth and effective radius in regions of the plume that are opaque to thermal IR measurements (AHI bands centred near 0.51, 0.64, 0.86, 1.6, 2.3 and 3.9 µm channels would be suitable; Prata and Grant2001; McGarragh et al.2015), retrieving other multi-layer scenarios (e.g. ice above ash or ash above ash) and improving the ash detection flag using machine learning (Picchiani et al.2011; Gray and Bennartz2015; Piontek et al.2021a).

Figure 13Ratios of volume extinction coefficients for different channel combinations as a function of effective radius. Volume extinction coefficients for 10.4, 11.2 and 12.4 µm are symbolised as k10, k11 and k12. The dashed lines with open circles indicate where the retrieval of the effective radius becomes ambiguous.


The ORAC retrievals provided here could either be used as a validation dataset for dispersion model simulations or incorporated into data assimilation schemes that require uncertainties at the grid/pixel level (e.g. Mingari et al.2022). Incorporation of these retrievals into such schemes could be used to develop quantitative now-casting/forecasting products that will aid VAACs in providing advice to airlines about quantitative ash concentrations in the future. In terms of implementing the retrieval scheme operationally, we have shown that one could use default a priori pressure settings in the troposphere and stratosphere and then select the run configuration with the lowest cost, on a per pixel basis, to identify whether or not ash is present in the stratosphere. Height estimates could then be refined as information from independent sources becomes available (e.g. from lidar, geostationary parallax and side-view heights).

Appendix A: Error propagation in the Mastin et al. (2009) equation

The relationship between the volumetric flow rate and the plume height has been developed over many years, as follows (Morton et al.1956; Wilson et al.1978; Settle1978; Bursik et al.1992; Sparks et al.1997; Mastin et al.2009):

(A1) H = a V ˙ b ,

where H is the plume height above vent level (km), V˙ is the volumetric flow rate (m3 s−1), and a and b are free parameters determined from an empirical fit. Equation (A1) can be rearranged to solve for the mass eruption rate, M˙ (kg s−1), as follows:

(A2) M ˙ = ρ d H a 1 / b ,

where ρd is the dense rock equivalent density of tephra (all particle sizes erupted from the volcano), which is typically assumed to be 2500 kg m−3 (Mastin et al.2009; Dioguardi et al.2020). Note that this density is distinct from the very fine ash particle density, ρ, which we assumed was 2300 ± 300 kg m−3 in the present study. If we assume all variables in Eq. (A2) are independent, then we can compute the uncertainty in M˙ from the partial derivatives as follows:

(A3) M ˙ ρ d = H a 1 / b ,

(A4) M ˙ H = 1 b H ρ d H a 1 / b ,

(A5) M ˙ a = - 1 b a ρ d H a 1 / b ,

(A6) M ˙ b = - ln ( H / a ) b 2 ρ d H a 1 / b .

The uncertainty in M˙ is then

(A7) ( α M ˙ ) 2 = M ˙ ρ d 2 ( α ρ d ) 2 + M ˙ H 2 ( α H ) 2 + M ˙ a 2 ( α a ) 2 + M ˙ b 2 ( α b ) 2 ,

which rearranges to the following:

(A8) α M ˙ M ˙ 2 = α ρ d ρ d 2 + 1 b α H H 2 + 1 b α a a 2 + ln ( H / a ) b α b b 2 = α ρ d ρ d 2 + 1 b 2 α H H 2 + α a a 2 + ln 2 ( H / a ) α b b 2 .

(A9) α M ˙ = ± M ˙ α ρ d ρ d 2 + 1 b 2 α H H 2 + α a a 2 + ln 2 ( H / a ) α b b 2 .

From Eq. (A8), we can see that the relative uncertainty in ρd is proportional to the relative uncertainty in M˙, whereas an increase in the relative uncertainty in H or a will result in a factor 1/b increase in the relative uncertainty in M˙. In addition, as the relative uncertainty in b increases, the relative uncertainty in M˙ will be an exponential function of H due to the ln(H/a)/b term and will thus dominate the uncertainty budget for large H (e.g. > 5.5 km for b=0.241).

To compute the total mass, MT, and its associated uncertainty, αMT, we must integrate M˙ with respect to time, t, as follows:

(A10) M T = t 0 t n M ˙ ( t ) d t .

The above integral can be approximated as a sum of mass eruption rates multiplied by discrete time periods (i.e. the temporal resolution of the satellite), Δt, so that, in the following:

(A11) M T i = 0 n M i ˙ Δ t = i = 0 n M i ,

where n is the total number of observation times, Mi˙ is the mass eruption rate at time, i, and Mi is the total mass at time, i, assuming constant Mi˙ over Δt (note that for AHI, Δt = 600 s). The uncertainty in the total mass can therefore be written as follows:


where we assume errors across each time step are uncorrelated. To evaluate M˙ and αM˙ at all observation times using Eqs. (A2) and (A9), we use the ORAC height time series shown in Fig. 5. Note that H is determined from the ORAC-retrieved height, hc, via H=hc-hv, where hv is the vent height of the Raikoke volcano (0.551 km a.s.l.). For ρd, a and b, we use the values provided in Mastin et al. (2009), which are 2500 kg m−3, 2.00 and 0.241, respectively. Uncertainty in H is taken from the ORAC retrievals. Uncertainty in ρd is not provided in Mastin et al. (2009), and so we conservatively assumed a relative uncertainty in ρd of 50 %. This value may seem high but, when propagated through Eqs. (A9) and (A13), the difference between a relative uncertainty of 10 % and 50 % converts to a difference in absolute uncertainty in the total mass of  1 Tg. Uncertainty associated with the a and b parameters is related to the sample size, eruption type and errors in plume height and volumetric flow rate taken from the cases used in the empirical fit. Based on the range of values reported and discussed in the literature (e.g. Bursik et al.1992; Sparks et al.1997; Mastin et al.2009), and noting the theoretical finding that M˙ is related to H to the fourth power (i.e. b=0.25Morton et al.1956), we assumed relative uncertainties in a and b of 90 % and 20 %, respectively.

Finally, the distal fine ash mass fraction, mf, is computed as follows:

(A14) m f = m T M T ,

where mT is the maximum total fine ash mass retrieved determined from ORAC. The associated uncertainty in mf is as follows:

(A15) α m f = ± m f α m T m T 2 + α M T M T 2 ,

where αmT is the uncertainty associated with the maximum total mass of very fine ash determined from ORAC.

Code and data availability

The ORAC code is open source and available from GitHub (, last access: 24 May 2022; Thomas et al.2009; Poulsen et al.2012; McGarragh et al.2018). The ORAC 10 min time series of plume-top heights and 1 h retrievals (all state variables and ash flag) out to 7 d after the eruption are available upon request from the corresponding author.


The supplement related to this article is available online at:

Author contributions

ATP, RGG and IAT developed the ideas for the paper. ATP led the writing of the paper, ran the satellite retrievals and conducted the data analysis. RGG, ACP, SRP, ATP and CP are developers of the ORAC software and made improvements to the code for the purposes of this study. All authors contributed to data interpretation and writing of the paper.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Special issue statement

This article is part of the special issue “Satellite observations, in situ measurements and model simulations of the 2019 Raikoke eruption (ACP/AMT/GMD inter-journal SI)”. It is not associated with a conference.


Roy G. Grainger and Isabelle A. Taylor have been supported by the NERC Centre for Observation and Modelling of Earthquakes, Volcanoes, and Tectonics (COMET). This study has been funded as part of NERC's support of the National Centre for Earth Observation (contract no. PR140015). Andrew T. Prata and Roy G. Grainger acknowledge funding from the NERC R4AsH project (grant no. NE/S003843/1). Isabelle A. Taylor and Roy G. Grainger acknowledge funding from the NERC V-Plus project (grant no. NE/S004025/1). We thank the two anonymous reviewers, for their detailed and thoughtful reviews, both of which helped to significantly improve the paper.

Financial support

This research has been supported by the Natural Environment Research Council (grant nos. NE/S003843/1 and NE/S004025/1) and the National Centre for Earth Observation (grant no. PR140015).

Review statement

This paper was edited by Christian von Savigny and reviewed by two anonymous referees.


Barber, C. B., Dobkin, D. P., and Huhdanpaa, H.: The quickhull algorithm for convex hulls, ACM Trans. Mathe. Softw., 22, 469–483,, 1996. a

Bessho, K., Date, K., Hayashi, M., Ikeda, A., Imai, T., Inoue, H., Kumagai, Y., Miyakawa, T., Murata, H., Ohno, T., Okuyama, A., Oyama, R., Sasaki, Y., Shimazu, Y., Shimoji, K., Sumida, Y., Suzuki, M., Taniguchi, H., Tsuchiyama, H., Uesawa, D., Yokota, H., and Yoshida, R.: An Introduction to Himawari-8/9 Japan's New-Generation Geostationary Meteorological Satellites, J. Meteorol. Soc. JPN Ser. II, 94, 151–183,, 2016. a

Bruckert, J., Hoshyaripour, G. A., Horváth, Á., Muser, L. O., Prata, F. J., Hoose, C., and Vogel, B.: Online treatment of eruption dynamics improves the volcanic ash and SO2 dispersion forecast: case of the 2019 Raikoke eruption, Atmos. Chem. Phys., 22, 3535–3552,, 2022. a

Bursik, M. I., Sparks, R. S. J., Gilbert, J. S., and Carey, S. N.: Sedimentation of tephra by volcanic plumes: I. Theory and its comparison with a study of the Fogo A plinian deposit, Sao Miguel (Azores), Bull. Volcanol., 54, 329–344,, 1992. a, b

Clarisse, L., Hurtmans, D., Prata, A. J., Karagulian, F., Clerbaux, C., De Maziére, M., and Coheur, P.-F.: Retrieving radius, concentration, optical depth, and mass of different types of aerosols from high-resolution infrared nadir spectra, Appl. Opt., 49, 3713,, 2010. a

Cooper, S. J., L'Ecuyer, T. S., Gabriel, P., Baran, A. J., and Stephens, G. L.: Objective Assessment of the Information Content of Visible and Infrared Radiance Measurements for Cloud Microphysical Property Retrievals over the Global Oceans. Part II: Ice Clouds, J. Appl. Meteorol. Climatol., 45, 42–62,, 2006. a

Corradini, S., Spinetti, C., Carboni, E., Tirelli, C., Buongiorno, M., Pugnaghi, S., and Gangale, G.: Mt. Etna tropospheric ash retrieval and sensitivity analysis using moderate resolution imaging spectroradiometer measurements, J. Appl. Remote Sens., 2, 023550,, 2008. a, b, c, d, e

Corradini, S., Montopoli, M., Guerrieri, L., Ricci, M., Scollo, S., Merucci, L., Marzano, F., Pugnaghi, S., Prestifilippo, M., Ventress, L., Grainger, R., Carboni, E., Vulpiani, G., and Coltelli, M.: A Multi-Sensor Approach for Volcanic Ash Cloud Retrieval and Eruption Characterization: The 23 November 2013 Etna Lava Fountain, Remote Sens., 8, 58,, 2016. a

Dacre, H. F., Grant, A. L. M., and Johnson, B. T.: Aircraft observations and model simulations of concentration and particle size distribution in the Eyjafjallajökull volcanic ash cloud, Atmos. Chem. Phys., 13, 1277–1291,, 2013. a

de Leeuw, J., Schmidt, A., Witham, C. S., Theys, N., Taylor, I. A., Grainger, R. G., Pope, R. J., Haywood, J., Osborne, M., and Kristiansen, N. I.: The 2019 Raikoke volcanic eruption – Part 1: Dispersion model simulations and satellite retrievals of volcanic sulfur dioxide, Atmos. Chem. Phys., 21, 10851–10879,, 2021. a

Deguine, A., Petitprez, D., Clarisse, L., Gudmundsson, S., Outes, V., Villarosa, G., and Herbin, H.: Complex refractive index of volcanic ash aerosol in the infrared, visible, and ultraviolet, Appl. Opt., 59, 884,, 2020. a

Dioguardi, F., Beckett, F., Dürig, T., and Stevenson, J. A.: The Impact of Eruption Source Parameter Uncertainties on Ash Dispersion Forecasts During Explosive Volcanic Eruptions, J. Geophys. Res.-Atmos., 125, e2020JD032717,, 2020. a

Folch, A., Mingari, L., Gutierrez, N., Hanzich, M., Macedonio, G., and Costa, A.: FALL3D-8.0: a computational model for atmospheric transport and deposition of particles, aerosols and radionuclides – Part 1: Model physics and numerics, Geosci. Model Dev., 13, 1431–1458,, 2020. a

Folch, A., Mingari, L., and Prata, A. T.: Ensemble-Based Forecast of Volcanic Clouds Using FALL3D-8.1, Front. Earth Sci., 9, 741841,, 2022. a

Francis, P. N., Cooke, M. C., and Saunders, R. W.: Retrieval of physical properties of volcanic ash using Meteosat: A case study from the 2010 Eyjafjallajökull eruption, J. Geophys. Res.-Atmos., 117, D00U09,, 2012. a, b, c, d, e, f, g, h, i

Gouhier, M., Eychenne, J., Azzaoui, N., Guillin, A., Deslandes, M., Poret, M., Costa, A., and Husson, P.: Low efficiency of large volcanic eruptions in transporting very fine ash into the atmosphere, Sci. Rep., 9, 1449,, 2019. a, b

Grainger, R. G., Peters, D. M., Thomas, G. E., Smith, A. J. A., Siddans, R., Carboni, E., and Dudhia, A.: Measuring volcanic plume and ash properties from space, Geological Society, London, Special Pub., 380, 293–320,, 2013. a

Gray, T. M. and Bennartz, R.: Automatic volcanic ash detection from MODIS observations using a back-propagation neural network, Atmos. Meas. Tech., 8, 5089–5097,, 2015. a

Gu, Y., Rose, W. I., and Bluth, G. J. S.: Retrieval of mass and sizes of particles in sandstorms using two MODIS IR bands: A case study of April 7, 2001 sandstorm in China, Geophys. Res. Lett., 30, 1805,, 2003. a

Gu, Y., Rose, W. I., Schneider, D. J., Bluth, G. J. S., and Watson, I. M.: Advantageous GOES IR results for ash mapping at high latitudes: Cleveland eruptions 2001, Geophys. Res. Lett., 32, L02305,, 2005. a

Harvey, N. J., Dacre, H. F., Webster, H. N., Taylor, I. A., Khanal, S., Grainger, R. G., and Cooke, M. C.: The Impact of Ensemble Meteorology on Inverse Modeling Estimates of Volcano Emissions and Ash Dispersion Forecasts: Grímsvötn 2011, Atmosphere, 11, 1022,, 2020. a

Harvey, N. J., Dacre, H. F., Saint, C., Prata, A. T., Webster, H. N., and Grainger, R. G.: Quantifying the impact of meteorological uncertainty on emission estimates and the risk to aviation using source inversion for the Raikoke 2019 eruption, Atmos. Chem. Phys., 22, 8529–8545,, 2022. a, b

Heidinger, A. K., Pavolonis, M. J., Holz, R. E., Baum, B. A., and Berthier, S.: Using CALIPSO to explore the sensitivity to cirrus height in the infrared observations from NPOESS/VIIRS and GOES-R/ABI, J. Geophys. Res., 115, D00H20,, 2010. a

Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz-Sabater, J., Nicolas, J., Peubey, C., Radu, R., Schepers, D., Simmons, A., Soci, C., Abdalla, S., Abellan, X., Balsamo, G., Bechtold, P., Biavati, G., Bidlot, J., Bonavita, M., Chiara, G., Dahlgren, P., Dee, D., Diamantakis, M., Dragani, R., Flemming, J., Forbes, R., Fuentes, M., Geer, A., Haimberger, L., Healy, S., Hogan, R. J., Hólm, E., Janisková, M., Keeley, S., Laloyaux, P., Lopez, P., Lupu, C., Radnoti, G., Rosnay, P., Rozum, I., Vamborg, F., Villaume, S., and Thépaut, J.: The ERA5 global reanalysis, Q. J. Roy. Meteorol. Soc., 146, 1999–2049,, 2020. a

Horváth, Á., Carr, J. L., Girina, O. A., Wu, D. L., Bril, A. A., Mazurov, A. A., Melnikov, D. V., Hoshyaripour, G. A., and Buehler, S. A.: Geometric estimation of volcanic eruption column height from GOES-R near-limb imagery – Part 1: Methodology, Atmos. Chem. Phys., 21, 12189–12206,, 2021a. a, b

Horváth, Á., Girina, O. A., Carr, J. L., Wu, D. L., Bril, A. A., Mazurov, A. A., Melnikov, D. V., Hoshyaripour, G. A., and Buehler, S. A.: Geometric estimation of volcanic eruption column height from GOES-R near-limb imagery – Part 2: Case studies, Atmos. Chem. Phys., 21, 12207–12226,, 2021b. a, b, c, d, e, f

Hyman, D. M. and Pavolonis, M. J.: Probabilistic retrieval of volcanic SO2 layer height and cumulative mass loading using the Cross-track Infrared Sounder (CrIS),, 2020. a

ICAO: Roadmap for International Airways Volcano Watch (IAVW) in Support of International Air Navigation, MOGVA ReferenceDocuments/IAVW Roadmap.pdf (last access: 29 November 2021), 2019. a

Inoue, T.: On the Temperature and Effective Emissivity Determination of Semi-Transparent Cirrus Clouds by Bi-Spectral Measurements in the 10 – Window Region, J. Meteorol. Soc. JPN, Ser. II, 63, 88–99,, 1985. a

Iwabuchi, H. and Hayasaka, T.: Effects of Cloud Horizontal Inhomogeneity on the Optical Thickness Retrieved from Moderate-Resolution Satellite Data, J. Atmos. Sci., 59, 2227–2242,<2227:EOCHIO>2.0.CO;2, 2002. a

Key, J. R.: Retrieval of cloud optical depth and particle effective radius at high latitudes using visible and thermal satellite data, p. 318, Paris, France,, 1995. a

Kylling, A., Kahnert, M., Lindqvist, H., and Nousiainen, T.: Volcanic ash infrared signature: porous non-spherical ash particle shapes compared to homogeneous spherical ash particles, Atmos. Meas. Tech., 7, 919–929,, 2014. a, b, c, d

Kylling, A., Kristiansen, N., Stohl, A., Buras-Schnell, R., Emde, C., and Gasteiger, J.: A model sensitivity study of the impact of clouds on satellite detection and retrieval of volcanic ash, Atmos. Meas. Tech., 8, 1935–1949,, 2015. a, b, c

Levenberg, K.: A method for the solution of certain non-linear problems in least squares, Q. Appl. Mathe., 2, 164–168,, 1944. a

Lindsey, D., Schmit, T. J., MacKenzie, W. M., Jewitt, C. P., Gunshor, M. M., and Grasso, L.: 10.35: atmospheric window on the GOES-R Advanced Baseline Imager with less moisture attenuation, J. Appl. Remote Sens., 6, 1,, 2012. a

Lu, S., Lin, H. X., Heemink, A., Segers, A., and Fu, G.: Estimation of volcanic ash emissions through assimilating satellite data and ground-based observations, J. Geophys. Res.-Atmos., 121, 10971–10994,, 2016. a

Mackie, S., Millington, S., and Watson, I. M.: How assumed composition affects the interpretation of satellite observations of volcanic ash: How assumed composition affects interpretation of ash observations, Meteorol. Appl., 21, 20–29,, 2014. a

Marquardt, D. W.: An Algorithm for Least-Squares Estimation of Nonlinear Parameters, J. Soc. Indust. Appl. Mathe., 11, 431–441,, 1963. a

Mastin, L., Guffanti, M., Servranckx, R., Webley, P., Barsotti, S., Dean, K., Durant, A., Ewert, J., Neri, A., Rose, W., Schneider, D., Siebert, L., Stunder, B., Swanson, G., Tupper, A., Volentik, A., and Waythomas, C.: A multidisciplinary effort to assign realistic source parameters to models of volcanic ash-cloud transport and dispersion during eruptions, J. Volcanol. Geothermal Res., 186, 10–21,, 2009. a, b, c, d, e, f, g, h, i, j, k, l

McGarragh, G. R., Thomas, G. E., Povey, A. C., Poulsen, C. A., and Grainger, R. G.: Volcanic ash retrievals using ORAC and satellite measurements in the visible and IR, Proc. ˜ATMOS 2015, Advances in Atmospheric Science and Applications, Heraklion, Greece, 8–12 June 2015 (ESA SP-735, November 2015), p. 8, 2015. a

McGarragh, G. R., Poulsen, C. A., Thomas, G. E., Povey, A. C., Sus, O., Stapelberg, S., Schlundt, C., Proud, S., Christensen, M. W., Stengel, M., Hollmann, R., and Grainger, R. G.: The Community Cloud retrieval for CLimate (CC4CL) – Part 2: The optimal estimation approach, Atmos. Meas. Tech., 11, 3397–3431,, 2018. a, b, c, d, e, f, g

McKee, K., Smith, C. M., Reath, K., Snee, E., Maher, S., Matoza, R. S., Carn, S., Mastin, L., Anderson, K., Damby, D., Roman, D. C., Degterev, A., Rybin, A., Chibisova, M., Assink, J. D., de Negri Leiva, R., and Perttu, A.: Evaluating the state-of-the-art in remote volcanic eruption characterization Part I: Raikoke volcano, Kuril Islands, J. Volcanol. Geotherm. Res., 419, 107354,, 2021. a, b, c, d, e, f, g, h

Mingari, L., Folch, A., Prata, A. T., Pardini, F., Macedonio, G., and Costa, A.: Data assimilation of volcanic aerosol observations using FALL3D+PDAF, Atmos. Chem. Phys., 22, 1773–1792,, 2022. a, b

Morton, B. R., Taylor, G. T., and Turner, J. S.: Turbulent gravitational convection from maintained and instantaneous sources, P. Roy. Soc. London Series A, 234, 1–23,, 1956. a, b

Muser, L. O., Hoshyaripour, G. A., Bruckert, J., Horváth, Á., Malinina, E., Wallis, S., Prata, F. J., Rozanov, A., von Savigny, C., Vogel, H., and Vogel, B.: Particle aging and aerosol–radiation interaction affect volcanic plume dispersion: evidence from the Raikoke 2019 eruption, Atmos. Chem. Phys., 20, 15015–15036,, 2020. a, b, c

Newman, S. M., Clarisse, L., Hurtmans, D., Marenco, F., Johnson, B., Turnbull, K., Havemann, S., Baran, A. J., and Haywood, J.: A case study of observations of volcanic ash from the Eyjafjallajókull eruption: 2. Airborne and satellite radiative measurements, p. 19,, 2012. a

Osborne, M. J., de Leeuw, J., Witham, C., Schmidt, A., Beckett, F., Kristiansen, N., Buxmann, J., Saint, C., Welton, E. J., Fochesatto, J., Gomes, A. R., Bundke, U., Petzold, A., Marenco, F., and Haywood, J.: The 2019 Raikoke volcanic eruption – Part 2: Particle-phase dispersion and concurrent wildfire smoke emissions, Atmos. Chem. Phys., 22, 2975–2997,, 2022. a, b, c, d, e, f

Pardini, F., Corradini, S., Costa, A., Esposti Ongaro, T., Merucci, L., Neri, A., Stelitano, D., and de Michieli Vitturi, M.: Ensemble-Based Data Assimilation of Volcanic Ash Clouds from Satellite Observations: Application to the 24 December 2018 Mt. Etna Explosive Eruption, Atmosphere, 11, 359,, 2020. a

Parol, F., Buriez, J. C., Brogniez, G., and Fouquart, Y.: Information Content of AVHRR Channels 4 and 5 with Respect to the Effective Radius of Cirrus Cloud Particles, J. Appl. Meteorol., 30, 973–984,, 1991. a

Pavolonis, M. and Sieglaff, J.: GOES-R Advanced Baseline Imager (ABI) Algorithm Theoretical Basis Document For Volcanic Ash (Detection and Height), p. 71, (last access: 17 October 2022), 2012. a

Pavolonis, M. J., Heidinger, A. K., and Sieglaff, J.: Automated retrievals of volcanic ash and dust cloud properties from upwelling infrared measurements: Retrieval of ash/dust cloud properties, J. Geophys. Res.-Atmos., 118, 1436–1458,, 2013. a, b, c, d, e, f, g, h, i, j, k, l, m

Pavolonis, M. J., Sieglaff, J., and Cintineo, J.: Spectrally Enhanced Cloud Objects - A generalized framework for automated detection of volcanic ash and dust clouds using passive satellite measurements: 2. Cloud object analysis and global application, J. Geophys. Res.-Atmos., 120, 7842–7870,, 2015a. a

Pavolonis, M. J., Sieglaff, J., and Cintineo, J.: Spectrally Enhanced Cloud Objects-A generalized framework for automated detection of volcanic ash and dust clouds using passive satellite measurements: 1. Multispectral analysis, J. Geophys. Res.-Atmos., 120, 7813–7841,, 2015b. a

Pavolonis, M. J., Sieglaff, J. M., and Cintineo, J. L.: Remote Sensing of Volcanic Ash with the GOES-R Series, in: The GOES-R Series, 103–124, Elsevier,, 2020. a

Picchiani, M., Chini, M., Corradini, S., Merucci, L., Sellitto, P., Del Frate, F., and Stramondo, S.: Volcanic ash detection and retrievals using MODIS data by means of neural networks, Atmos. Meas. Tech., 4, 2619–2631,, 2011. a

Piontek, D., Bugliaro, L., Kar, J., Schumann, U., Marenco, F., Plu, M., and Voigt, C.: The New Volcanic Ash Satellite Retrieval VACOS Using MSG/SEVIRI and Artificial Neural Networks: 2. Validation, Remote Sens., 13, 3128,, 2021a. a

Piontek, D., Hornby, A., Voigt, C., Bugliaro, L., and Gasteiger, J.: Determination of complex refractive indices and optical properties of volcanic ashes in the thermal infrared based on generic petrological compositions, J. Volcanol. Geotherm. Res., 411, 107174,, 2021b. a

Pouget, S., Bursik, M., Webley, P., Dehn, J., and Pavolonis, M.: Estimation of eruption source parameters from umbrella cloud or downwind plume growth rate, J. Volcanol. Geotherm. Res., 258, 100–112,, 2013. a

Poulsen, C. A., Siddans, R., Thomas, G. E., Sayer, A. M., Grainger, R. G., Campmany, E., Dean, S. M., Arnold, C., and Watts, P. D.: Cloud retrievals from satellite data using optimal estimation: evaluation and application to ATSR, Atmos. Meas. Tech., 5, 1889–1910,, 2012. a, b, c, d, e, f

Prabhakara, C., Fraser, R. S., Dalu, G., Wu, M.-L. C., Curran, R. J., and Styles, T.: Thin Cirrus Clouds: Seasonal Distribution over Oceans Deduced from Nimbus-4 IRIS, J. Appl. Meteorol., 27, 379–399,<0379:TCCSDO>2.0.CO;2, 1988. a

Prata, A. J.: Infrared radiative transfer calculations for volcanic ash clouds, Geophys. Res. Lett., 16, 1293–1296,, 1989a. a, b, c

Prata, A. J.: Observations of volcanic ash clouds in the 10–12 µm window using AVHRR/2 data, Int. J. Remote Sens., 10, 751–761,, 1989b. a

Prata, A. J. and Grant, I. F.: Retrieval of microphysical and morphological properties of volcanic ash plumes from satellite data: Application to Mt Ruapehu, New Zealand, Q. J. Roy. Meteorol. Soc., 127, 2153–2179,, 2001. a, b, c, d, e, f

Prata, A. J. and Prata, A. T.: Eyjafjallajókull volcanic ash concentrations determined using Spin Enhanced Visible and Infrared Imager measurements, J. Geophys. Res.-Atmos., 117, D00U23,, 2012. a, b, c, d, e, f

Prata, A. T., Siems, S. T., and Manton, M. J.: Quantification of volcanic cloud top heights and thicknesses using A-train observations for the 2008 Chaitén eruption, J. Geophys. Res.-Atmos., 120, 2928–2950,, 2015. a

Prata, A. T., Young, S. A., Siems, S. T., and Manton, M. J.: Lidar ratios of stratospheric volcanic ash and sulfate aerosols retrieved from CALIOP measurements, Atmos. Chem. Phys., 17, 8599–8618,, 2017a. a, b

Prata, A. T., Folch, A., Prata, A. J., Biondi, R., Brenot, H., Cimarelli, C., Corradini, S., Lapierre, J., and Costa, A.: Anak Krakatau triggers volcanic freezer in the upper troposphere, Sci. Rep., 10,, 2020. a

Prata, A. T., Mingari, L., Folch, A., Macedonio, G., and Costa, A.: FALL3D-8.0: a computational model for atmospheric transport and deposition of particles, aerosols and radionuclides – Part 2: Model validation, Geosci. Model Dev., 14, 409–436,, 2021. a, b, c, d

Prata, F., Woodhouse, M., Huppert, H. E., Prata, A., Thordarson, T., and Carn, S.: Atmospheric processes affecting the separation of volcanic ash and SO2 in volcanic eruptions: inferences from the May 2011 Grímsvötn eruption, Atmos. Chem. Phys., 17, 10709–10732,, 2017b. a, b

Prata, G. S., Ventress, L. J., Carboni, E., Mather, T. A., Grainger, R. G., and Pyle, D. M.: A New Parameterization of Volcanic Ash Complex Refractive Index Based on NBO/T and SiO 2 Content, J. Geophys. Res.-Atmos., 124, 1779–1797,, 2019. a, b, c

Reed, B. E., Peters, D. M., McPheat, R., and Grainger, R. G.: The Complex Refractive Index of Volcanic Ash Aerosol Retrieved From Spectral Mass Extinction, J. Geophys. Res.-Atmos., 123, 1339–1350,, 2018. a, b

Rodgers, C. D.: Retrieval of atmospheric temperature and composition from remote measurements of thermal radiation, Rev. Geophys., 14, 609,, 1976. a

Rodgers, C. D.: Inverse methods for atmospheric sounding: Theory and Practice, Series on Atmospheric, Oceanic and Planetary Physics, Vol. 2, World Scientific, Singapore, 240 pp., ISBN: 9789810227401, 2000. a, b

Saunders, R., Hocking, J., Turner, E., Rayer, P., Rundle, D., Brunel, P., Vidot, J., Roquet, P., Matricardi, M., Geer, A., Bormann, N., and Lupu, C.: An update on the RTTOV fast radiative transfer model (currently at version 12), Geosci. Model Dev., 11, 2717–2737,, 2018. a

Schneider, D. J., Rose, W. I., Coke, L. R., Bluth, G. J. S., Sprod, I. E., and Krueger, A. J.: Early evolution of a stratospheric volcanic eruption cloud as observed with TOMS and AVHRR, J. Geophys. Res.-Atmos., 104, 4037–4050,, 1999. a

Sears, T. M., Thomas, G. E., Carboni, E., A. Smith, A. J., and Grainger, R. G.: SO 2 as a possible proxy for volcanic ash in aviation hazard avoidance, J. Geophys. Res.-Atmos., 118, 5698–5709,, 2013. a

Settle, M.: Volcanic eruption clouds and the thermal power output of explosive eruptions, J. Volcanol. Geotherm. Res., 3, 309–324,, 1978. a

Smirnov, S., Nizametdinov, I., Timina, T., Kotov, A., Sekisova, V., Kuzmin, D., Kalacheva, E., Rashidov, V., Rybin, A., Lavrenchuk, A., Degterev, A., Maksimovich, I., and Abersteiner, A.: High explosivity of the June 21, 2019 eruption of Raikoke volcano (Central Kuril Islands); mineralogical and petrological constraints on the pyroclastic materials, J. Volcanol. Geotherm. Res., 418, 107346,, 2021. a, b

Soda, R.: Infrared Absorption Spectra of Quartz and Some other Silica Modification, B. Chem. Soc. JPN, 34, 1491–1495,, 1961. a

Sparks, R. S. J., Bursik, M. I., Carey, S. N., Gilbert, J., Glaze, L. S., Sigurdsson, H., and Woods, A. W. (Eds.): Volcanic Plumes, Wiley, 574 pp., ISBN: 9780471939016, 1997. a, b

Stamnes, K., Tsay, S.-C., Wiscombe, W., and Laszlo, I.: DISORT, a general-purpose Fortran program for discrete-ordinate-method radiative transfer in scattering and emitting layered media: documentation of methodology, Tech. rep., Dept. of Physics and Engineering Physics, Stevens Institute of Technology Hoboken, NJ 07030, 112 pp., (last access: 17 October 2022), 2000. a

Stohl, A., Prata, A. J., Eckhardt, S., Clarisse, L., Durant, A., Henne, S., Kristiansen, N. I., Minikin, A., Schumann, U., Seibert, P., Stebel, K., Thomas, H. E., Thorsteinsson, T., Tørseth, K., and Weinzierl, B.: Determination of time- and height-resolved volcanic ash emissions and their use for quantitative ash dispersion modeling: the 2010 Eyjafjallajökull eruption, Atmos. Chem. Phys., 11, 4333–4351,, 2011. a

Thomas, G. E., Poulsen, C. A., Sayer, A. M., Marsh, S. H., Dean, S. M., Carboni, E., Siddans, R., Grainger, R. G., and Lawrence, B. N.: The GRAPE aerosol retrieval algorithm, Atmos. Meas. Tech., 2, 679–701,, 2009. a, b

Van Eaton, A. R., Amigo, Ã., Bertin, D., Mastin, L. G., Giacosa, R. E., González, J., Valderrama, O., Fontijn, K., and Behnke, S. A.: Volcanic lightning and plume behavior reveal evolving hazards during the April 2015 eruption of Calbuco volcano, Chile, Geophys. Res. Lett., 43, 3563–3571,, 2016. a

Vicente, G. A., Davenport, J. C., and Scofield, R. A.: The role of orographic and parallax corrections on real time high resolution satellite rainfall rate distribution, Int. J. Remote Sens., 23, 221–230,, 2002. a

Wang, C., Yang, P., Baum, B. A., Platnick, S., Heidinger, A. K., Hu, Y., and Holz, R. E.: Retrieval of Ice Cloud Optical Thickness and Effective Particle Size Using a Fast Infrared Radiative Transfer Model, J. Appl. Meteorol. Climatol., 50, 2283–2297,, 2011. a

Watts, P., Mutlow, C., Baran, A., and Zavody, A.: Study on cloud properties derived from Meteosat Second Generation observations, EUMETSAT ITT, 97, 181, (last access: 17 October 2022), 1998. a, b, c

Watts, P. D., Bennartz, R., and Fell, F.: Retrieval of two-layer cloud properties from multispectral observations using optimal estimation, J. Geophys. Res., 116, D16203,, 2011. a

Wen, S. and Rose, W. I.: Retrieval of sizes and total masses of particles in volcanic clouds using AVHRR bands 4 and 5, J. Geophys. Res., 99, 5421,, 1994. a, b, c, d, e, f, g

Western, L. M., Watson, M. I., and Francis, P. N.: Uncertainty in two-channel infrared remote sensing retrievals of a well-characterised volcanic ash cloud, B. Volcanol., 77, 67,, 2015. a

Wilkins, K. L., Mackie, S., Watson, M., Webster, H. N., Thomson, D. J., and Dacre, H. F.: Data insertion in volcanic ash cloud forecasting, Ann. Geophys., 57, 1–6,, 2015. a

Wilkins, K. L., Watson, I. M., Kristiansen, N. I., Webster, H. N., Thomson, D. J., Dacre, H. F., and Prata, A. J.: Using data insertion with the NAME model to simulate the 8 May 2010 Eyjafjallajókull volcanic ash cloud, J. Geophys. Res.-Atmos., 121, 306–323,, 2016. a, b

Wilson, L., Sparks, R. S. J., Huang, T. C., and Watkins, N. D.: The control of volcanic column heights by eruption energetics and dynamics, J. Geophys. Res.-Solid Earth, 83, 1829–1836,, 1978. a

Winker, D. M., Vaughan, M. A., Omar, A., Hu, Y., Powell, K. A., Liu, Z., Hunt, W. H., and Young, S. A.: Overview of the CALIPSO Mission and CALIOP Data Processing Algorithms, J. Atmos. Ocean. Technol., 26, 2310–2323,, 2009. a

Winker, D. M., Liu, Z., Omar, A., Tackett, J., and Fairlie, D.: CALIOP observations of the transport of ash from the Eyjafjallajókull volcano in April 2010, J. Geophys. Res.-Atmos., 117, D00U15,, 2012. a

Witham, C., Hort, M., Thomson, D., Devenish, B., Webster, H., and Beckett, F.: The current volcanic ash modelling set- up at the London VAAC, p. 11, (last access: 17 October 2022), 2019.  a

Yamanouchi, T., Suzuki, K., and Kawaguchi, S.: Detection of Clouds in Antarctica from Infrared Multispectral Data of AVHRR, J. Meteorol. Soc. JPN Ser. II, 65, 949–962,, 1987. a

Yang, P., Feng, Q., Hong, G., Kattawar, G. W., Wiscombe, W. J., Mishchenko, M. I., Dubovik, O., Laszlo, I., and Sokolik, I. N.: Modeling of the scattering and radiative properties of nonspherical dust-like aerosols, J. Aerosol Sci., 38, 995–1014,, 2007. a

Yu, T., Rose, W. I., and Prata, A. J.: Atmospheric correction for satellite-based volcanic ash mapping and retrievals using “split window” IR data from GOES and AVHRR, J. Geophys. Res., 107, D16,, 2002. a

Zidikheri, M. J. and Lucas, C.: A Computationally Efficient Ensemble Filtering Scheme for Quantitative Volcanic Ash Forecasts, J. Geophys. Res.-Atmos., 126, e2020JD033094,, 2021. a

Short summary
Satellite observations are often used to track ash clouds and estimate their height, particle sizes and mass; however, satellite-based techniques are always associated with some uncertainty. We describe advances in a satellite-based technique that is used to estimate ash cloud properties for the June 2019 Raikoke (Russia) eruption. Our results are significant because ash warning centres increasingly require uncertainty information to correctly interpret, aggregate and utilise the data.