Optimal ash particle refractive index model for simulating the brightness temperature spectrum of volcanic ash clouds from satellite infrared sounder measurements

Using data from the Infrared Atmospheric Sounding Interferometer (IASI) measurements of volcanic ash clouds and radiative transfer calculations, we identify the optimal refractive index model for simulating the measured brightness temperature spectrum of volcanic ash material. We assume that the optimal refractive index model has the smallest root mean square of the brightness temperature difference between measurements and simulations for channels in the wavenumber range 10 of 750–1400 cm−1 and compare 21 refractive index models for optical properties of ash particles, including recently published models. From the results of numerical simulations for 164 pixels of IASI measurements for ash clouds from 11 volcanoes, we found that the measured brightness temperature spectrum could be well simulated using certain newly established refractive index models. In the cases of Eyjafjallajökull and Grímsvötn ash clouds, the optimal refractive index models determined through numerical simulation correspond to those deduced from the chemical composition of ash samples for the same volcanic 15 eruption events. This finding suggests that infrared sounder measurement of volcanic ash clouds is an effective approach to estimating the optimal refractive index model. However, discrepancies between the estimated refractive index models based on satellite measurements and the associated volcanic rock types were observed for some volcanic events.


Introduction
For continuous near-real time monitoring of volcanic ash clouds (VACs) during both day and night, satellite observations 20 obtained using infrared sensors are particularly useful (Mackie and Watson, 2014;Mackie et al., 2016;Clarisse and Prata, 2016). Because volcanic silicate has a characteristic absorption band around 10 μm due to the stretching vibration of the Si-O bond, VACs containing fine ash particles can be easily detected and distinguished from water or ice clouds through measurement of the brightness temperature difference between two channels in this absorption band, such as 10.8 μm and 12 μm (Wen and Rose, 1994;Prata and Grant, 2001). Applying further constraints related to the atmospheric profile and 25 microphysical properties of the ash particles, such as the particle size distribution, the physical properties of VACs (cloud height, particle effective radius, optical depth, and the associated ash mass loading) were determined using measurements in other infrared channels from high-resolution satellite imagers (Francis et al., 2012;Pavolonis et al., 2013Pavolonis et al., , 2015aPavolonis et al., , 2015bPrata and Lynch, 2019) and satellite infrared sounders Ventress et al., 2016;Clarisse and Prata, 2016).
These retrieval methods are based on the results of radiative transfer calculations for the measured infrared brightness 30 temperature of VACs under assumed atmospheric profiles and surface conditions, therefore the optical properties of the ash particles are essential. In particular, the complex refractive index (RI) is an important upstream parameter for estimation of ash optical properties and their wavelength dependence. A typical RI model over a wide range of infrared wavelengths, the "andesite" model proposed by Pollack et al. (1973) (hereinafter PL1973) has been adopted for radiative transfer calculations for volcanic ash analysis. Although the andesite 35 model of PL1973 provides suitable approximations for analysis of some VACs (Pavolonis et al., 2013), simulations of the infrared brightness temperature spectrum (BTS) of VACs based on radiative transfer calculations from PL1973 andesite often caused large discrepancies relative to calculations based on satellite infrared sounder measurements (Gangale et al., 2010; https://doi.org/10.5194/amt-2021-103 Preprint. Discussion started: 20 April 2021 c Author(s) 2021. CC BY 4.0 License.
Atmospheric Radiative Transfer Simulator (Buehler et al., 2010). Upward radiance of a sounder channel at TOA is calculated as the weighted sum of number of monochromatic radiances .
where is the weighting coefficient �1 = ∑ �. Using a diverse atmospheric profile dataset to compare the calculated with values derived from exact integration of channel response functions using line-by-line calculations, the monochromatic 85 wavenumber and weight were determined, as was the number , using a simulated annealing technique. To simulate the radiance of the 8461 IASI channels, a total of 14,132 monochromatic wavenumber values were determined. The rootmean-square error of for clear-sky conditions was less than the IASI noise equivalent brightness temperature difference, except in the channels of the O3 band at ≈ 1050 cm −1 and the CO2 band at 2200 ≤ ≤ 2400 cm −1 .
A plane-parallel single homogeneous layer VAC was assumed for calculations of multiple scattering by cloud particles. The 90 values of and in Eq. (1) determined under clear-sky conditions can be applied to calculations in a cloudy atmosphere with high accuracy (Holl et al., 2012), and we used the same formulation of Eq. (1) for radiance calculations in cloudy conditions caused by volcanic ash. For monochromatic radiative transfer calculations at infrared wavelengths, we originally formulated an analytical expression of TOA radiance in which the terms for cloud multiple scattering were described by coefficients depending only on cloud geometry and cloud optical properties. A look-up-table (LUT) of coefficients for various ash cloud 95 conditions was prepared for calculation of monochromatic radiance in advance using the discrete ordinate method with 64 streams. In our benchmark test, the CPU time for TOA radiance calculation using MBCRM over 100 layers of atmosphere and 8461 IASI channels was 0.73 (0.34) sec under cloudy (clear) conditions.

Ash cloud parameters for radiative transfer calculations and RI models 100
Using LUTs for ash optical properties, numerical simulations of IASI measurements of VACs were conducted using MBCRM.
A prolate spheroid shape with a 1:2 axis ratio and a log-normal size distribution (Hansen and Travis, 1974) with a geometric standard deviation of 0.74 (Pavolonis et al., 2013) were assumed for ash particles. The optical properties of the particles were calculated using the T-Matrix method (Mishchenko et al., 2002) for particle effective radii from 0.1 to 20 μm (particle size is defined as the radius of a sphere with equivalent volume). In satellite infrared sounder measurements, volcanic ash 105 mainly affects TOA brightness temperature in the wavenumber range of ≤ 1400cm −1 , and LUTs of volcanic ash optical properties were prepared for 640 cm −1 ≤ ≤ 1400 cm −1 at an interval of ∆ = 10 cm −1 . The optical properties of each IASI channel were estimated through linear interpolation of the LUT. TOA radiance and brightness temperature were calculated under the assumption that the ash clouds are in thermal equilibrium with the ambient air.
For RI models of volcanic ash particles, we used the datasets of Reed et al. (2018) and by Prata et al. (2019) (hereinafter 110 RE2018 and PG2019). Additionally, we prepared datasets using the two RI models MP_A and MP_R. These models were originally described as "andesite" and "rhyolite (obsidian)" models by PL1973, and artificial weak absorption features were added to them. The RI models applied in this work are listed in Table 1, and the characteristics of the models are described below. Using the assumptions of the Rayleigh continuous distribution of ellipsoids (CDE) scattering model for the sampled ash particles and the Lorentz formulation of the RI to ensure consistency with the Kramers-Kronig relationship, Reed et al. (2018) provided a spectral RI dataset of typical volcanic ash samples at wavelengths of 0.33-19 µm (RE2018). RI data for eight ash samples from seven volcanoes (Askja, Aso, Eyjafjallajökull-(a) and -(b), Grímsvötn, Nisyros, Spurr, and Tongariro) are available from the Aerosol Refractive Index Archive (http://eodg.atm.ox.ac.uk/ARIA/index.html) as well as in the supporting 150 information of PG2019 (https://doi.org/10.1029/2018JD028679). The spectral RIs of RE2018 over the wavenumber range of 650 cm −1 ≤ ≤ 1400 cm −1 are plotted in Fig. 1. In the wavenumber range of 700 cm −1 ≤ ≤ 800 cm −1 , the RI models of RE2018 showed weak absorption features that are not clearly expressed in the andesite and rhyolite models of PL1973. In previous attempts to simulate sounder-observed BTSs, the calculated BTSs tended to have positive bias in this wavenumber range (Gangale et al., 2010;Clarisse et al., 2013;Ishimoto et al., 2016;Clarisse and Prata, 2016). With the RE2018 RI models, 155 improvements of BTS fitting are expected, although the relationships of chemical composition and volcano type (mafic or felsic) with the absorption index (imaginary part of the complex refractive index: ) at 700 cm −1 ≤ ≤ 800 cm −1 remain unclear. All RI models of the RE2018 dataset have relatively weak absorption peaks ( ≤ 1.0) and similar wavenumber dependence of between 1110 and 1250 cm −1 . The individual RI models of RE2018, RE020 (Aso) and RE080 (Tongariro), show similar wavenumber dependence of , although their NBO/T and SiO2 wt. % values differ significantly. In addition, 160 RE050 (Grímsvötn) has a broad feature of that differs from other RI models of RE2018.  Table 1

RI models parameterised by Prata et al. (2019): PG2019
Using the RE2018 RI dataset and the results of laboratory analysis of the chemical composition of the ash samples, Prata et al.
(2019) provided a new parameterisation for the ash RI based on SiO2 content and NBO/T. According to , 185 NBO/T and SiO2 wt. % in the bulk composition of the ash samples show stronger correlations with the derived RI over a broad spectral range compared to glass composition. Furthermore, parameterisation using NBO/T provides better results than SiO2 wt. % in both visible and infrared wavelengths. In this study, we used 11 RI models parameterised using NBO/T between 0 and 1 with a step size of 0.1. The ( , ) data for corresponding NBO/T values were calculated using a supporting spreadsheet of , which is available from the publisher's website. Compared to the RE2018 dataset, RI models using 190 NBO/T parameterisation show simplified wavelength dependences in both the real and imaginary parts of the RI due to the use of linear regression. However, this simple parameterisation is useful for VAC retrieval based on satellite remote sensing because an appropriate RI model can be specified in advance based on the assumed chemical composition of the ash material.
In that context, it is essential to verify how well PG2019 RI models reproduce the BTSs of infrared sounder measurements.

Modified RI models based on andesite and rhyolite of Pollack et al. (1973): MP_A and MP_R
The two RI models MP_A and MP_R are based on the andesite and rhyolite (obsidian from Little Glass Mt. California) models of PL1973, respectively, and we artificially applied features of weak absorption in the wavenumber range of 700-850 cm −1 to the absorption index of each model, referring to the results of BTS analysis reported by (Gangale et al., 2010;Ishimoto et 210 al., 2016). These RI models were used for comparisons to the results of BTS simulations using RI models based on the RE2018/PG2019 datasets.

IASI ash cloud data 230
To identify the optimal RI models by comparing measured and simulated BTSs, continuous data from an infrared sounder that covers the wavenumber range of 650-1400 cm −1 is desirable, and VAC data collected by IASI were used in this work. The IASI Level-1C granule data were acquired from the National Oceanic and Atmospheric Administration Comprehensive Large Array-Data Stewardship System (https://www.avl.class.noaa.gov/saa/products/). As a primary condition for BTS analysis, the brightness temperature difference between 10.7 μm (934.5 cm −1 ) and 12.2 μm (819.5 cm −1 ) ( ≡ 10.7 − 12.2 ) 235 was used, and pixels with < −2 K for VACs over the ocean were selected. For discriminating the results of BTS simulations from different RI models using radiative transfer calculations, VAC pixels with more negative values of are better, as the BTSs of such pixels shows distinct wavenumber dependence due to absorption in the infrared window region.
The value −2 K was chosen as the minimum requirement for discriminating the difference in BTSs between the results of different RI models. We excluded VAC data over land due to large uncertainties in spectral surface emissivity and surface 240 temperature.
Our radiative transfer model MBCRM conducts forward calculations for a single layer of VAC with no contamination from meteorological clouds (MCs) consisting of water or ice particles. Because MC contamination in VAC pixels can cause large estimation errors for ash cloud parameters and RI models (Kylling et al., 2015), contaminated pixels must be excluded to ensure the validity of the estimated RI model based on the results of our BTS analysis. In the case of VAC covered with dense 245 MC layers, most measured pixels are rejected due to the condition < −2 K. On the other hand, if a portion of the pixel area is covered by MC or if the VAC is located above the MC layers, the MC can cause a negative bias in the measured BTS and still meet < −2 K. As a result, large errors may occur for the estimation of ash cloud parameters under the assumption of a single VAC layer. Furthermore, the detailed wavelength dependence of the BTS for such contaminated pixels may differ from that of the pure VAC. For these reasons, we used VAC data from daytime IASI measurements, and pixels that 250 were presumed to be affected by MC were excluded through reference to visible true-colour imagery for the same area from the Moderate Resolution Imaging Spectroradiometer (MODIS) onboard the Terra and Aqua satellites. MODIS images were provided by the Level-1 and Atmosphere Archive & Distribution System Distributed Active Archive Center, and are available online (https://atmosphere-imager.gsfc.nasa.gov/images/l1b-granules). The visible true-colour image product of the Japanese geostationary meteorological satellite Himawari-8 Advanced Himawari Imager (AHI) (Bessho et al., 2016) was also used for 255 volcanic eruptions since 2015 in the coverage area of Himawari-8 (images were downloaded from the Japan Aerospace Exploration Agency Himawari Monitor homepage: https://www.eorc.jaxa.jp/ptree/index.html). Because we used daytime IASI measurements that could be referenced to MODIS images, the number of VAC pixels available for our analysis was significantly reduced. Furthermore, the evaluation of IASI pixels using MODIS visible images does not always exclude MC contamination due to the time difference between IASI and MODIS measurements as well as the difference in sensor spatial 260 resolution. Therefore, we added another condition for IASI pixels in our analysis. As the temperature of the MC layer is generally lower than that of the sea surface at the same geographic location, a VAC above an MC layer tends to have a lower infrared brightness temperature than a VAC with no MC if the cloud parameters of the VAC are the same. Using brightness temperature at wavenumber , where the measured brightness temperature is maximised in the region of 750 cm −1 ≤ ≤ 900 cm −1 due to a local minimum of the absorption index , the difference between the brightness temperatures of 265 clear sky ( ) and ash clouds ( ) was defined as [ The value of was derived from VAC measurements and ( ) was calculated using our forward model. Then, a threshold of was set according to the value of . Figure 4 shows a plot of ( ) with respect to in 1171 IASI pixels of VAC measurements for the volcanic eruptions listed in Table 2. The pixels are selected based on the primary condition < −2K over ocean areas, including night-time data. The clear-sky brightness temperature ( ) was calculated using 270 atmospheric profiles, sea surface temperature, and atmospheric pressure at the surface from the results of global analysis (GANAL) via numerical modelling by the Japan Meteorological Agency. In Fig. 4, pixels are discriminated between Eyjafjallajökull (blue) and other volcanic eruptions (red), and the two datasets in the plot show similar distribution patterns. If a pixel of pure VAC observation is contaminated with MC, the negative value of is likely to become smaller, while the value of ∆ is likely to increase. Therefore, the upper right data in Fig. 4 may have higher probabilities of MC 275 contamination than the lower left data on the plot. In our retrieval analysis, a threshold value of ∆ (K) was defined to reduce the fraction of MC-contaminated pixels: and ∆ < ∆ . ℎ was applied as an additional retrieval condition. The coefficients of Eq.
(2) were artificially determined, and approximately 35% of all pixels with < −2K were rejected due to this condition for the 280 Eyjafjallajökull VAC shown in Fig. 4. The aim of this condition is to reduce the fraction of MC-contaminated pixels in our analysis; this threshold cannot be used to determine whether an individual pixel is contaminated with MC. Even for a pure VAC of a single homogeneous layer, ∆ may be greater than ∆ . ℎ depending on the cloud parameters of the VAC. In addition to MC contamination, the negative value of decreases and ∆ increases for VACs with larger particle sizes and optical depths, as shown in the scatter plot of Fig. 4, which is essentially the same as the split-window plotting 285 method (Prata, 1989;Wen and Rose, 1994;Prata and Grant, 2001)  For retrieval analysis, we selected up to 10 VAC pixels with large negative values of from a granular dataset of IASI measurements. Radiative transfer calculations and BTS simulations for the selected VAC pixels were conducted using 310 assumptions for three VAC parameters, pressure height of cloud top , optical depth at wavelength 11 μm ( = 907.488 cm −1 ), and effective radius of ash particles , which were retrieval variables, and using the temperature/water vapour profiles and sea surface temperature/pressure from GANAL as fixed values. We set the pressure height of the cloud base to = + 100 hP a for simplicity. For the ozone profile, the total column of ozone was treated as a retrieval variable, and the relative values for the initial ozone profile from GANAL were fixed. A homogeneous volcanic SO2 gas layer 315 with a pressure difference of 100 hP a between the top and base was assumed, with the pressure height of the layer top and total column used as variables. The estimated values of the retrieval variables were those for which the root mean square (RMS) between measured and calculated brightness temperatures ( and ) in the assumed channels of number N is minimised. The optimal RI models were generally identified by comparing the RMS values of various RI models.
(3) 320 Retrieval calculations were conducted using the scattering properties for ash particles derived from the selected RI model and data for IASI channels of 650 cm −1 ≤ ≤ 1400 cm −1 with the following procedure: a) VAC parameters � , , � are estimated using and for channels excluding the ozone band of 980 cm −1 ≤ ≤ 1080 cm −1 and SO2 bands of 1100 cm −1 ≤ ≤ 1210 cm −1 and 1320 cm −1 ≤ ≤ 1395 cm −1 .
b) The total columns of ozone and SO2 are estimated using the channels 980 cm −1 ≤ ≤ 1080 cm −1 for ozone and 325 1320 cm −1 ≤ ≤ 1395 cm −1 for SO2, assuming the VAC parameters � , , � as fixed values. The wavenumbers 1100 cm −1 ≤ ≤ 1210 cm −1 are excluded from estimation of the SO2 column because the brightness temperature of these channels is strongly influenced by the applied ash RI model.
c) The RMS for the channels 750 cm −1 ≤ ≤ 1400 cm −1 is then calculated, considering the estimated VAC, ozone, and SO2 parameters as fixed. We excluded the range of 650 cm −1 ≤ < 750 cm −1 to avoid RMS values related to error in the 330 GANAL atmospheric profiles. The date and number of pixels in the selected IASI data for eruptions from 11 volcanoes are listed in Table 2. The average value of the wavenumber ̅ at the local maximum of brightness temperature and SiO2 wt. % of tephra samples from the literature are included for reference. In particular, is related to the minimum of the absorption index , and has been reported 335 as an important parameter for identifying ash material from infrared sounder measurements (Gangale et al., 2010;Clarisse and Prata, 2016). Using a total of 21 RI models, which are listed in Table 1, the minimum RMS of brightness temperature in the wavenumber range 750 cm −1 ≤ ≤ 1400 cm −1 was determined for each RI model. In the following sections, we discuss the results for the best fitting (optimal) RI models determined from RMS values for the volcanic events in Table 2. The detailed results for minimal RMS and VAC parameters � , , �, ash top height ℎ estimated from 340 and the GANAL profile, and total column SO2 for each IASI pixel are provided in the Supplemental Material. In this study, we assume that the ash material ejected from a given volcano remains the same throughout the period of each eruption event in Table 2, and also assume that the optimal RI model is the one for which the sum of the RMS for all pixels of a given volcano is smallest. Even if the total RMS was the smallest, we rejected an RI model as unrealistic if the converged values of were the smaller than 0.2 μm for many pixels in the retrieval analyses. In addition to the differences among RI models, 345 the size of ash particles effectively changes the wavelength dependence of the simulated BTS (Clarisse and Prata, 2016;Ishimoto et al., 2016), and small values of tend to be estimated when a relatively mafic RI model is used with low SiO2 wt. % or large NBO/T values. The relationship between the estimated value and the RI model in the retrieval analysis is discussed in section 5.1.
350 Table 2: Volcanic eruption events used for BTS analysis of IASI measurements and the results of the optimal RI models. The optimal RI models are identified through comparison between measured and simulated BTSs for IASI pixels on the date "Date." Two RI models in the column "Optimal RI model" indicate that the results of BTS fitting were almost the same. The average wavenumber

Refractive index models estimated through IASI measurement simulation
Selected IASI measurements for the eruption events of 11 volcanoes and the results of the optimal RI models are listed in 375 Table 2. Ranking plots for the results of RMS analysis for each measurement pixel are shown in Fig. 5. Ranking based on the smallest RMS was conducted for the 21 RI models. Detailed results for each RI model in each pixel, including the minimum RMS, estimated ash cloud parameters � , ℎ , , �, and SO2 column, are provided in the Supplemental Material. The RI models from the PG2019 dataset were selected as optimal for 10 volcanic events. Furthermore, high-ranking RI models tended to converge at a specific value of NBO/T in the PG2019 dataset. This result suggests that NBO/T, which can be 380 determined from the chemical composition of the ash, is an effective indicator of the appropriate RI model for satellite VAC analysis. On the other hand, the andesite model of PL1973 (MP_A) was not selected as the best RI model despite having a better RMS than the original andesite model of PL1973. This finding clearly shows that high-accuracy BTS simulations can be achieved using the novel RI models of RE2018 and PG2019. The results of VAC analysis for individual volcanoes are discussed in the following subsections. Typical results of the BTS simulations are shown in the figures, and the estimated VAC 385 parameters, SO2 column, and RI model used for each simulation are summarised in Table 3.    Fig. 6. For these 460 IASI measurements, true-colour images from the MODIS product ( Fig. 6a and d) near the IASI measurement time were obtained. Through comparison with the colour image shown in Fig. 6a, the influence of MCs appears small in the vicinity of IASI pixels with highly negative values in Fig. 6b. The results of BTS simulation with RE030 and PG040 were very similar and agreed well with the measured BTS, although some differences occurred for the estimated VAC parameters �ℎ , , � (Table 3). Because PG040 is parameterised by NBO/T = 0.4 and this NBO/T value is close to that of the 465 Eyjafjallajökull-(a) ash sample at 0.38, the wavenumbers of the local minimum and local maximum of the absorption index obtained using PG040 are similar to those of RE030 (Figs. 1 and 2). Thus, the results of BTS simulations shown in Fig. 6c arise from the similarity of the spectral pattern of between PG040 and RE030. Another Eyjafjallajökull VAC was located more than 1000 km from the volcano on 9 May 2010 ( Fig. 6d-f) Table 3.

510
Using the retrieval results for 50 pixels of Eyjafjallajökull VAC measurements under various RI models (PG070, PG040, and PG010), the dependence of the RI model on the estimated VAC parameters �ℎ , , � was investigated (Fig. 7a-c).
Although the estimated VAC parameters differed among the RI models tested, no systematic bias was apparent for ℎ or .
On the other hand, the results for show a strong dependence on RI model selection (Fig. 7b). Within the PG2019 dataset, 515 small values of tend to be estimated for mafic RI models with large NBO/T. By contrast, relatively large values were derived for felsic RI models. Similar dependence of RI models on estimated was also confirmed for our VAC retrievals related to other volcanic events listed in Table 2. Using a pixel with small differences in retrieved ℎ and among the three RI models (pixel no. 30), the results of BTS simulations with different values using PG070 and PG010 are plotted in Fig.   7d and 7e, respectively. If we assume that the RI model of PG040 is adequate for Eyjafjallajökull ash material and that the 520 retrieval result = 1.07 μm using PG040 is close to the true value, a positive bias of the simulated BTS in the region of 900 cm −1 ≤ ≤ 1230 cm −1 occurs, except within the ozone band, when we use a more mafic RI model (PG070) than PG040 with the same value (blue line of Fig. 7d). The positive bias of the BTS is mitigated by applying a smaller value of in the BTS simulation (green line). By contrast, a negative bias of the BTS occurs in the same wavenumber range (blue line of Fig. 7e) when we use a more felsic RI model (PG010), and fitting of the calculated BTS improves when a larger value is 525 used for the ash particles (green line of Fig. 7e). Thus, the retrieval results of are sensitive to the wavenumber dependence of the particle absorption property around 900 cm −1 ≤ ≤ 1230 cm −1 . Therefore, selection of the proper RI model is essential, especially for the estimation of . This result also suggests that RI model selection may have a strong influence on the estimation of ash column density, as it is directly related to the retrieved (Pavolonis et al., 2013). In our BTS analysis of the Eyjafjallajökull VAC, PG040 and RE030 were selected as the optimal RI models as their fitting results were better than 530 those of other RI models in wavenumber regions other than 900 cm −1 ≤ ≤ 1230 cm −1 (in particular, 750 cm −1 ≤ ≤ 900 cm −1 ).

Grimsvötn
In addition to Eyjafjallajökull, an RI model from the laboratory experiments of RE2018 and IASI measurements of ash clouds were obtained for the Grímsvötn eruptions in May 2011. From IASI data for the Grímsvötn VAC, one pixel at 12:08 UTC on 580 22 May and 10 pixels at 11:47 UTC on 23 May were selected for our BTS analysis. The measured IASI pixels of < −2K on 23 May (Fig. 8b) were estimated to have low MC contamination based on the MODIS colour image taken at 12:05 UTC (Fig. 8a). Moreover, from this MODIS colour image, MCs in the area of the VAC were determined to be located at higher altitudes than the VAC (Prata et al., 2017b), suggesting that the MC-contaminated pixels were effectively excluded from the target pixels in our BTS analysis. 585 https://doi.org/10.5194/amt-2021-103 Preprint. Discussion started: 20 April 2021 c Author(s) 2021. CC BY 4.0 License.
For the VAC of Grímsvötn eruptions, the BTS based on IASI measurements was reported to produce good simulation results using the basalt RI model of PL1973 (Newman et al., 2012). The results of our BTS simulations showed that a small RMS was obtained for the PG070 RI model using the PG2019 dataset, and RI models derived from Grímsvötn (RE050) and Spurr (RE070) ash samples within RE2018 dataset returned relatively small RMS values. Absorption features related to SO2 gas were not present in the BTSs of the 11 analysed pixels, in accordance with the report of SO2 gas separation from the ash clouds 590 (Prata et al., 2017b), and thus PG070 and RE050 were selected as the optimal RI models based on their total RMS scores. The estimate of NBO/T for the Grímsvötn ash sample is 0.74-0.75, which is the highest value within the RE2018 dataset . Because the PG2019 dataset is derived from NBO/T parameterisation using the RE2018 dataset, the wavelength dependence of PG070 is similar to that of RE050. As a result, similar VAC parameters were estimated for PG070 and RE050, and the average values over 11 pixels �ℎ � , ̅ , ̅ � = (1.58 km, 0.93 μm, 2.06) were almost identical for these two RI 595 models. This result for generally agrees with published results for a Grímsvötn ash sample (1.1 μm) (Reed et al., 2018). The best RI models for simulating the observed BTSs of the Eyjafjallajökull and Grímsvötn VACs had similar NBO/T values 620 to ash samples from the same volcanic eruptions. This result provides evidence that the BTS of volcanic ash clouds obtained from satellite infrared sounder measurements is closely related to the spectrum of the ash RI and cloud parameters. Furthermore, it supports the proposals of Reed et al. (2018) and  that an appropriate RI model for VAC retrieval from satellites can be determined using the NBO/T (or SiO2 content) of the ash samples. In addition, our results suggest that the optimal RI model can be identified using only infrared sounder measurements for the VAC under certain atmospheric 625 conditions to conduct radiative transfer calculations, provided that enough RI models have been prepared in advance. https://doi.org/10.5194/amt-2021-103 Preprint. Discussion started: 20 April 2021 c Author(s) 2021. CC BY 4.0 License.

Calbuco
In the series of Calbuco volcanic activity from 22 April 2015, IASI data for brightness temperature at 13:29-13:32 UTC on 630 24 April were used for our retrieval analysis. At this time, ash clouds with negative were observed over the Atlantic Ocean near the estuary of the La Plata River at latitudes of 35-38° S and longitudes of 55-56° W (Calbuco_A: Fig. 9b), as well as over the Pacific Ocean near the west coast of Chile at 32-33° S and 73-76° W (Calbuco_B: Fig. 9e). From the MODIS colour image taken at 13:35 UTC (Fig. 9a), it appeared that the effect of MCs on the BTS was small for pixels with large negative values of in Calbuco_A. Although another VAC was confirmed around latitude 40° S in the MODIS image 635 shown in Fig. 9a, the IASI pixels had slightly negative or positive values, suggesting ice particles within or above the VAC layer. BTS analyses were conducted for 15 pixels with large negative values in Calbuco_A. Small total RMS values were obtained from the PG010-PG030 RI models with the PG2019 dataset and RE040 with the RE2018 dataset, and the smallest RMS was obtained from RE040 (Fig. 5). The estimated values were small overall, and were less than 0.2 μm for most pixels when relatively mafic RI models (PG030-PG100 of PG2019 dataset) were applied. For the RE040 RI model, 640 low of less than 0.2 μm were obtained for most of the analysed pixels, and the smallest effective radius = 0.1 μm in our LUT was estimated for 10 of the 15 total pixels within Calbuco_A. These small values are mainly related to the fitting of the BTS in the range of 1070 cm −1 ≤ ≤ 1230 cm −1 , as discussed in section 5.1 (Fig. 7). The area of pixels with large negative in Fig. 9b is ~1500 km from Calbuco volcano, and the particle size of the VAC should be reduced due to depositional segregation during transportation in the atmosphere. Although a < 0.2 μm may not be too small for the ash 645 particles present in this area, we considered the PG020 RI model, which had the second smallest RMS after RE040 and averaged ≥ 0.2 μm, to be the optimal RI model in this study. As shown in Fig. 9c, the absorption feature of SO2 gas is clearly apparent in the measured BTS, and the results of BTS simulations approximately agree with the measurements.
However, the RMS values between measurements and simulations were relatively large for all pixels and systematic deviation from the measurements occurred even with PG020. 650 MODIS observed the Calbuco_B region at 15:00 UTC (Fig. 9d). We inferred the presence of an optically thin VAC in an area with little MCs, although the MODIS observation was made 1.5 hours after that of IASI. As SO2 absorption features were not present in the BTS from IASI measurements, we deduced that the measured BTS patterns mainly arise from the optical properties of the VAC. In BTS simulations for 10 pixels in Calbuco_B, high-ranking RI models tended to have greater NBO/T in the PG2019 dataset compared with the results of Calbuco_A (Fig. 5); results from the RI models PG020-PG080 and RE040-655 RE050 had small RMS values and the smallest RMS was obtained for PG050. The reason for this difference between the results of Calbuco_A and Calbuco_B has not been determined. According to the compositional analyses of Romero et al. (2016) and Deguine et al. (2020), Calbuco ash samples are of the basaltic andesite type with 55-56.3 SiO2 wt. %, which is a similar or more mafic composition than Eyjafjallajökull ash samples. Our results for Calbuco_B pixels are relatively consistent with the compositional analysis results, but the PG050 RI had large RMS errors for Calbuco_A pixels 660 in our retrieval analysis. By contrast, the difference of RMS values between PG050 and PG020 is generally small in BTS simulations, as shown in Fig. 9f. Therefore, to explain the BTS based on both Calbuco_A and Calbuco_B measurements, we treated PG020 as the optimal RI model for Calbuco VAC retrieval.

Kirishimayama and Nishinoshima
During the Kirishimayama eruption on 26 January 2011, ash plumes were transported in the southeast direction over the Pacific Ocean, and IASI measured the VAC to the south of the Japanese main islands at 00:08 UTC on 27 January (Fig. 10b). The 700 MODIS image taken at 01:05 UTC (Fig. 10a) suggests that the VAC in the region of highly negative is distributed above the streaks of MCs, and VACs with lower MC contamination are expected in locations between the MC streaks. Due to the 1-hour time difference between IASI and MODIS measurements, the pixels for analysis were determined using only IASI data. Among 30 pixels of < −2K, 21 pixels were rejected under the condition of Eq.
(2), suggesting that many IASI pixels were contaminated with MCs. For the BTSs of the remaining nine pixels, the RI models RE080, RE020, and PG020 705 provided small RMS values. The results of BTS simulations for RE080 and RE020 generally agree with measured values (Fig.   10c). However, some small differences between the measured and calculated BTSs were confirmed, such as a negative bias in the RE080 simulation around 920 cm −1 and a positive bias in the RE020 simulation around 900-1000 cm −1 . From the ranking plot shown in Fig. 5 for the PG2019 dataset, the Kirishimayama VAC at this time tended to fit the results of RI models with smaller NBO/T values than the Eyjafjallajökull and Grímsvötn ash clouds. This result is reasonable considering the 710 reported SiO2 contents of the ash samples and the average local maximum of ̅ ≈ 830 cm −1 , which is a greater wavenumber https://doi.org/10.5194/amt-2021-103 Preprint. Discussion started: 20 April 2021 c Author(s) 2021. CC BY 4.0 License. than the local maxima for Eyjafjallajökull and Grímsvötn (Table 2). However, the possibility of MC contamination of the BTSs for all VAC pixels in this area cannot be excluded. Further precise analysis of the Kirishimayama VAC, including nighttime measurements with no MC contamination, is necessary. An algorithm for detecting MC contamination in the VAC using infrared channels will be essential to solving this problem. 715 Compared to Kirishimayama, the VACs of Nishinoshima, which were observed by IASI at 00:08 (Fig. 10e) and 23:47 UTC on 31 July 2020, were less affected by MC contamination due to the simultaneous measurements taken by the AHI onboard the geostationary satellite Himawari-8 (Fig. 10d). In BTS simulations for 18 pixels with large negative values of , the RI models RE070 with the RE2018 dataset and PG060-PG080 with the PG2019 dataset returned small RMS values and small effective radius estimates ( ≤ 0.6 μm) for all pixels. The smallest RMS was obtained with RE070. As shown in Fig. 10f, 720 the calculated BTS fit well to the measured BTS, especially at wavenumbers ≥ 1050 cm −1 , but systematic errors in the range of 750 cm −1 ≤ ≤ 1000 cm −1 were also observed. In their analysis, a mixture of these two RI models was assumed as the a priori RI and the spectral RI in the wavenumber range 760 of ν ≥ 1100 cm −1 was omitted due to missing spectral data between 1137 and 1216 cm −1 from the AIRS measurements.
Furthermore, MC contamination was not considered for pixels with highly negative values, although MC contamination may cause large errors in RI estimation. For the Kelut VAC at this time, Kylling (2016) reported that a mixture of ice clouds and ash clouds in an IASI pixel can explain the observed BTS from the PL1973 andesite model. In this work, retrieval analyses of the Kelut VAC were repeated using IASI pixels with less MC contamination, as estimated from visible 765 images and from the condition of Eq. (2).
In the IASI data of the Kelut VAC observed at 02:08 UTC on 14 February (Fig. 11c), all pixels with < −2K in the region of 8.0° S ≤ ≤ 9.5° S and 106.5° E ≤ ≤ 109.5° E were rejected under the condition of Eq.
(2), and four pixels (indicated by the circle in Fig. 11c) were retained for our retrieval analysis. This result suggests that most VAC measurements are affected by MC contamination, which is consistent with the conclusion of Kylling (2016) that IASI measures both ash and 770 ice clouds in the same pixel. The area of the four retained pixels is presumed to be covered with thin ash clouds and have less MCs based on the MODIS true-colour image taken at 03:35 UTC (Fig. 11a) and the Himawari-7 visible monochromatic image taken at 02:00 UTC (Fig. 11b). The results of BTS simulations showed that the BTS of the Kelut VAC could be well simulated using relatively felsic RI models with small NBO/T values and the PG2019 dataset (Fig. 5), with PG000 and PG010 leading to the smallest RMS values (Fig. 11d). Although the measured BTSs for pixels in the area of 8.0° S ≤ ≤ 9.5° S and 775 106.5° E ≤ ≤ 109.5° E are negatively biased due to MC contamination, their spectral features, such as the wavenumber of the local maximum ( ), were similar to those shown in Fig. 11d (plots of BTSs for the Kelut VAC are also shown in Clarisse and Prata (2016)). Whether these measured BTSs can be simulated well by internal or external mixing between the andesite model of VAC and ice cloud remains unclear. In addition to internal mixing between ash and ice based on effective medium theory, a simple external mixing process between individual ash and ice layers was examined through preliminary 780 calculations, and the BTS simulations were unsuccessful. We could not conclude that Kelut ash clouds were composed of felsic ash materials based on the chemical composition (~55 SiO2 wt. %) of Kelut ash samples from the same eruption (Maeno et al., 2019). However, the measured BTSs of the Kelut VACs in this eruption event can be well explained by assuming that ash particles in the Kelut VACs had a similar infrared RI feature to PG000 and PG010.

Puyehue-Cordón Caulle (PCC)
From the beginning of the eruption on 4 June 2011, large amounts of volcanic ash and SO2 gas were ejected from PCC, and the plumes reached an altitude of 12-13 km. Subsequently, ash clouds remained at high latitudes of the southern hemisphere for a long time (Theys et al., 2013;Klüser et al., 2013;Clarisse et al., 2013). BTS analysis was conducted for ash clouds observed over the Atlantic Ocean at 10:41-10:47 UTC and 12:29 UTC on 6 June. In the area of the IASI measurements at 825 10:41-10:47 UTC (Fig. 12c), MODIS observations were made at 11:05 UTC (Fig. 12b) and 12:45 UTC (Fig. 12a). From these MODIS images, ash clouds in the area with highly negative values of based on IASI measurements were distributed above broken MCs, and BTS simulations were performed for 29 IASI pixels that met the condition of Eq. (2). In addition, ten additional IASI pixels were selected in the area of latitude 40-42° S and longitude 60-63° W from the data obtained at 12:29 UTC to avoid low-level dense MCs (circles in Fig. 12e, f). As shown in the ranking plot in Fig. 5, almost all pixels have similar 830 spectral features, and smaller RMS values were obtained for RI models with smaller NBO/T values in the PG2019 dataset.
The smallest RMS was obtained for PG000 among all RI models. In accordance with the results of BTS analysis conducted by Newman et al. (2012), felsic RI models agreed well with the compositional analysis results of ash sampled by Castro et al. (2013). Comparing BTSs from simulations and measurements, the simulation results obtained using PG000 differ significantly from the observations (Fig. 12d, g). Although the simulation results show a good fit to measurements in the wavenumber 835 region of 820-1000 cm −1 and at wavenumbers greater than 1130 cm −1 , a large difference in brightness temperature at wavenumbers less than 820 cm −1 suggests an excess of relative absorption in the PG000 RI model. Moreover, the simulated BTS values were significantly larger than measurements between 1070 and 1130 cm −1 , particularly for pixels in the area https://doi.org/10.5194/amt-2021-103 Preprint. Discussion started: 20 April 2021 c Author(s) 2021. CC BY 4.0 License. shown in Fig. 12c. In the results of VAC parameter retrieval using PG000, the effective radius was less than 0.5 μm for most pixels, and tended to decrease as distance from PCC increased, and most of the retrieval results showed ≤ 840 0.3 μm for pixels around ( , ) = (27° S, 25° W) (circle in Fig. 12c), which is 4500-4600 km from PCC. Considering a decrease in large particles due to deposition during a long transport period in the atmosphere, ≤ 0.3 μm may be realistic in the case of PCC. Ash clouds from PCC eventually circled the southern hemisphere three times. Very small ash particles may increase the longevity of PCC ash clouds (Carn and Krotkov, 2016). According to the results of measurements by Cloud-Aerosol Lidar with Orthogonal Polarization onboard Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations, PCC 845 ash clouds are mainly composed of fine-mode non-spherical ash particles with a low sulphate contribution (Vernier et al., 2013;Prata et al., 2017a). For BTSs in the range of 1070-1200 cm −1 , the RI model MP_R provides a better fit than PG000 (data not shown). This result suggests that an RI model with a stronger absorbing feature in the index around 1050 cm −1 than that of PG000 is better for explaining the measured BTS of PCC ash clouds.

Bezymianny, Zhupanovsky, Rinjani, and Sarychev Peak
In analyses of VACs from Bezymianny, Zhupanovsky, Rinjani, and Sarychev Peak, a small number of IASI pixels could be obtained due to the relatively small areas of the VACs (Bezymianny, Zhupanovsky, Rinjani), or due to distinct MC 880 https://doi.org/10.5194/amt-2021-103 Preprint. Discussion started: 20 April 2021 c Author(s) 2021. CC BY 4.0 License. contamination of most < −2K pixels in MODIS visible images (Sarychev Peak). Therefore, we omit detailed discussion of these volcanic events, and simply present the results of our analyses. The dates and times of the satellite observations and retrieval parameters are listed in Tables 2 and 3, respectively. For the VAC of Sarychev Peak, one IASI pixel in an area of low MC contamination based on comparison with MODIS imagery (indicated by a circle in Fig. 13j) was used for retrieval analysis. For these four volcanic events, the results of our analyses indicate good BTS fits between simulations 885 and measurements, but whether the best-fit RI model is the optimal RI model remains unclear due to the small number of retrievals. Validation of the results using other eruption events is necessary in the future.

Summary and discussions
Using infrared brightness temperature data from IASI measurements of VACs and radiative transfer calculations, we aimed to identify optimal RI models for ash particles that could reproduce the measured BTSs. We applied a total of 21 RI models to 940 BTS simulations, including eight RI models using the RE2018 dataset, 11 RI models using the PG2019 dataset, and two additional RI models based on the andesite and rhyolite models of PL1973. For the eruption events of 11 volcanoes, 164 daytime IASI pixels containing apparent volcanic ash features ( ≤ −2K) with low MC contamination were selected with reference to MODIS and Himawari-8/AHI visible-colour images and an additional brightness temperature condition. The BTS simulation for each pixel included retrieval of VAC parameters �ℎ , , �, SO2 column content, and the RMS 945 between the measured and simulated BTSs for each RI model was calculated. The results of the RI models were ranked based on the smallest RMS values for each pixel, and the optimal RI models for the target volcanoes were determined on the basis of the total RMS. Typical BTS simulation results were discussed for the ash clouds of some volcanoes. The findings obtained from these analyses are summarised below.
Using RI models based on new RI datasets (RE2018, PG2019), the accuracy of brightness temperature simulation for VACs 950 in the atmospheric window region was greatly improved. In particular, RI models using the PG2019 dataset performed well, and they were selected as the optimal RI models for 10 out of 11 volcanic events. This finding suggests that the NBO/T parameterisation of  is effective for the RI at infrared window wavelengths. Moreover, in analyses of each pixel, RI models parameterised with NBO/T often showed better BTS fits than the RI models of Reed et al. (2018) , which were estimated directly from samples of specific volcanoes. This result indicates that simplified and generalised RI models 955 created by removing small fluctuation components provide a convenient method to simulate measured BTS features.
In our BTS simulations, RI models for Eyjafjallajökull (RE030, RE040) with the RE2018 dataset and RI models parameterised using similar NBO/T values to those of Eyjafjallajökull ash samples (PG030, PG040) represented the measured BTSs of Eyjafjallajökull ash clouds with good accuracy. Such consistency was also observed for Grímsvötn ash clouds. The one-toone correlation between ash samples and ash cloud measurements based on RI for the same volcanic event indicates that the 960 optimal RI model for ash materials can be determined from hyperspectral sounder measurements as well as NBO/T and SiO2 wt. % data obtained from compositional analysis of ash samples under the condition that sufficient RI models are available in advance. Although we used the RI models of the RE2018 and PG2019 datasets, which were developed from laboratory experiments of RI and the chemical compositions of ash samples, our basic approach is to estimate the RI model from satellite measurements, and the results for Eyjafjallajökull and Grímsvötn support the validity of our process. https://doi.org/10.5194/amt-2021-103 Preprint. Discussion started: 20 April 2021 c Author(s) 2021. CC BY 4.0 License.
The RI models developed from BTS simulations of IASI measurements were not always consistent with models deduced from the reported chemical compositions of ash samples. For ash clouds released from Calbuco, Kelut, Zhupanovsky, and Bezymianny, the RI models determined from IASI measurements were more felsic (or had lower NBO/T) than the rock type or SiO2 wt. % of tephra samples of the corresponding volcano. Furthermore, some RI models were excluded from the optimal RI models despite good BTS fitting results because the retrieved effective radii were smaller than 0.2μm. Negative values of 970 brightness temperature difference are greater for ash particles of smaller sizes, and therefore our analysis is biased toward ash clouds with small . In addition, the retrieved values tended to be smaller when mafic RI models were used compared to the results of felsic RI models. This difference is due to fitting of the measured BTS in the wavenumber range of 950 cm −1 ≤ ≤ 1230 cm −1 , excluding the ozone absorption band. An effective radius of less than 0.2 μm is too small for fresh ash clouds shortly after eruption, but may be realistic for ash clouds that have undergone long-distance transport in the 975 atmosphere. Moreover, the condensation of volcanic sulphate on pure ash particles during transport might alter the inherent optical properties of the ash clouds. In particular, a retrieved ≤ 0.2 μm may be reasonable for ash clouds measured far from the volcano, such as those from Calbuco and Puyehue-Cordón Caulle observed over the Atlantic Ocean.
We used a modified version of the andesite model (MP_A), in which a weak absorption feature in the range of 700-850 cm −1 was added. Although the modified RI model showed better fitting results than the original andesite model of PL1973 in our 980 BTS analysis, it was not selected as the best RI model for all IASI pixels in this study. This result indicates that substantial improvement in satellite retrieval results for VACs can be expected from replacing the conventional andesite model with a proper RI model based on the RE2018 and PG2019 datasets. We explored the reproducibility of measured BTSs over the entire wavenumber range of 650-1400 cm −1 , including bands for CO2, water vapour, ozone, and SO2, and noted that the estimated RI models may not produce good results for VAC parameters in retrievals from multi-channel satellite imagers using certain 985 infrared wavelength channels.
In this study, we used MODIS and Himawari-8 daytime visible images for rough estimation of MC contamination in IASI measurement areas. For that reason, the number of IASI pixels available for analysis was significantly reduced, and statistical evaluation of the retrieval results for some volcanic events was difficult. To achieve a sufficient number of analyses, thorough evaluation of MC contamination of night-time VAC data from infrared sounder measurements is necessary, and methods for 990 such evaluation should be explored in future research. In recent years, plans for mounting infrared sounders on geostationary satellites have been undertaken. High-frequency observations of distinct volcanic plumes increase the number of BTS data available, and precise retrieval is expected with further improvements in VAC analysis to combine infrared sounder and highresolution imager measurements.
We tested 21 RI models in this article, and conspicuous discrepancies in measured and calculated BTSs at specific wavenumber 995 ranges were confirmed in all RI models for some volcanic events. This finding suggests that the number of RI models is still insufficient for simulating all BTS patterns of VACs for various types of volcanoes. Recently, another RI dataset based on laboratory analysis of ash samples of six volcanoes was published by Deguine et al. (2020), and the data are available from the journal website. We conducted the same BTS analysis for VAC pixels for the Calbuco, Eyjafjallajökull, Grímsvötn, and Puyehue-Cordón Caulle events using the corresponding RI models. However, the results of these BTS simulations did not 1000 change the ranking or selection of optimal RI models in this paper. Compared to the RE2018 and PG2019 datasets, peaks of the absorption index in the RI models of Deguine et al. (2020) are shifted toward longer wavelengths, which appears to increase the RMS between measurements and our simulation results. Deguine et al. (2020) notes that this difference in the absorption index is related to the differing assumed shape of ash particles in laboratory analysis of ash samples, as Deguine et al. (2020) considered a spherical shape and applied Mie theory to the optical properties of sampled ash particles, while Reed 1005Reed et al. (2018 used the results of the Rayleigh CDE. If the spectral features of the resultant RI model are affected by the assumed ash shape, realistic shape models for ash particles may lead to novel RI models that agree well with the BTSs of ash clouds https://doi.org/10.5194/amt-2021-103 Preprint. Discussion started: 20 April 2021 c Author(s) 2021. CC BY 4.0 License.