Unified quantitative observation of coexisting volcanic sulfur dioxide and sulfate aerosols using ground-based Fourier transform infrared spectroscopy

We developed an optimal-estimation algorithm to simultaneously retrieve, for the first time, coexisting volcanic gaseous SO₂ and sulfate aerosols (SA) from ground-based Fourier transform infrared (FTIR) observations. These effluents, both linked to magmatic degassing process and subsequent atmospheric evolution processes, have overlapping spectral signatures leading to mutual potential interferences when retrieving one species without considering the other. We show that significant overestimations may be introduced in SO₂ retrievals if the radiative impact of coexistent SA is not accounted for, which may have impacted existing SO₂ long-term series, e.g. from satellite platforms. The method was applied to proximal observations at Masaya volcano, where SO₂ and SA concentrations, and SA acidity, were retrieved. A gas-to-particle sulfur partitioning of 400 and a strong SA acidity (sulfuric acid concentration: 65 %) were found, consistent with past in situ observations at this volcano. This method is easily exportable to other volcanoes to monitor magma extraction processes and the atmospheric sulfur cycle in the case of ash-free plumes.

1 Introduction

Volcanic gas and particulate emissions affect tropospheric and stratospheric compositions, air quality and the environment, the distribution and optical properties of low and high clouds, the Earth radiation budget from the regional to the global scale, and therefore climate (e.g. von Glasow et al., 2009; Robock, 2000). One of the most important environmental pollutants and the main source of radiative forcing from volcanoes is long-lived acidic and highly reflective sulfate aerosols (SA), directly emitted (primary SA) or formed by gas-to-particle conversion of sulfur dioxide (SO₂) emissions (secondary SA).

Observing these volcanic emissions and their atmospheric processes and variability in space and time using ground-based and satellite remote sensing is a crucial step towards understanding and quantifying their environmental and climatic impacts. Proximal integrated observations of different and interacting gaseous and particulate volcanic effluents is also important to gain insights into magmatic degassing processes and eruption forecasting. Fourier transform infrared (FTIR) spectroscopy (e.g. Oppenheimer et al., 1998; Francis et al., 1998; Duffell et al., 2003; Burton et al., 2007) is a powerful tool in this context. Ground-based FTIR is an ideal instrument for measuring magmatic degassing as it allows for remote quantification of the major magmatic gases, including H₂O, CO₂, SO₂, HCl, HF and also trace gases such as SiF₄, CO and OCS. Typically a radiation source is needed, and this may be an artificial infrared lamp or hot rock or lava (e.g. Allard et al., 2005). Correspondingly, open-path (OP) FTIR is based on an atmospheric path of a known distance between an artificial or natural source and an FTIR optical system. These spectrometers collect spectra over a large spectral range, which contains broad spectral features arising from volcanic aerosols as well as more finely structured molecular absorption signatures. However, until now, these broadband aerosol features have not been utilised as the research focus has been exclusively on magmatic gas emissions. This means a very rich resource of information on volcanic aerosols and their processes is potentially available in previously measured FTIR spectra from volcanoes worldwide.

The possibility of sporadic detection of SA, without specific physicochemical characterisation, from high-spectral-resolution infrared satellite instruments has been shown in the past for relatively strong volcanic eruptions with stratospheric injection (Haywood et al., 2010; Karagulian et al., 2010). In this work, we demonstrate that systematic detection, quantification and chemical characterisation of volcanic sulfate aerosols using OP-FTIR spectrometry is feasible. As SO₂ and SA emissions are inter-connected by magmatic degassing processes and by the subsequent short-term and small-scale atmospheric processes (e.g. Sellitto et al., 2017a) and have absorption features in the same spectral range (e.g. Sellitto and Legras, 2016; Sellitto et al., 2017b), we explore the possibility of characterising simultaneously these two volcanic effluents. This co-retrieval is intended to limit mutual biases, which are expected if SO₂ and SA are retrieved individually (e.g. Guermazi et al., 2017). In addition, this co-retrieval provides, in principle, simultaneous information on two interacting species, contributing a constraint to the inner volcanic and atmospheric sulfur cycle. This new inversion method (described in Sect. 3) is then applied to OP-FTIR observations taken at Masaya volcano during a measurement campaign conducted in 1998 (described in Sect. 2). Results are given and discussed in Sect. 4 and conclusions are drawn in Sect. 5.

2 The case study: Masaya volcano and the observation geometry

Masaya volcano (11.98^∘ N, 86.16^∘ W) is located in the Central American volcanic belt, which runs from Guatemala to Costa Rica in the north–south direction. It is situated about 25 km south-east of Managua city, Nicaragua, with an elevation of about 600 m. It is structured as a basaltic-andesitic shield caldera. Masaya is one of the world's most persistent sources of passive magmatic degassing with a relatively stable degassing rate and episodic stronger events (Rymer et al., 1998). The active vent releases SO₂, whose fluxes range from about 500 to 2500 t d⁻¹ (Duffell et al., 2003; Mather et al., 2006; de Moor et al., 2013; Carn et al., 2017). Masaya is, therefore, one of the largest contributors of volcanic gas emissions in the Central American arc (de Moor et al., 2017). One of Masaya's most remarkable features is its accessibility; it is literally a drive-in volcano with a car park next to the main degassing crater. This means that it has been used as a natural laboratory to conduct investigations into magmatic degassing and aerosol evolution. In situ air mass sampling and subsequent laboratory analyses of Masaya's plume to investigate aerosol composition and burden were conducted during several campaigns (e.g. Allen et al., 2002; Martin et al., 2011).

The data used in this work were collected not with direct sampling but with OP-FTIR remote sensing during a two-year (1998–1999) measurement campaign. In particular, we focus on high-quality spectra collected during a measurement session from about 16:40 to 17:10 UTC on 15 March 1998. During these measurements, Masaya was not in an eruptive period and observations are therefore representative of the normal passive degassing activity. During the campaign, a MIDAC Corporation OP-FTIR spectrometer was deployed on the top of Santiago crater. The observation geometry for this instrumental set-up is depicted in Fig. 2 of Horrocks et al. (1999): an infrared lamp and an FTIR spectrometer were placed on two sides of the active Masaya crater so that the radiation of the lamp is observed by the FTIR after it is transmitted through the plume with a total distance of 518 m. More details on the FTIR system and its technical specifications, the campaign, and its major results can be found in (Horrocks et al., 1999; Burton et al., 2000).

3 Methodology

With reference to Fig. 2 of Horrocks et al. (1999), the radiance spectrum measured by the spectrometer when the plume is in the line of sight is

$$\begin{array}{}\text{(1)}& I\left(\mathit{\lambda }\right)=I_{\mathrm{lamp}}\left(\mathit{\lambda }\right)e^{-{\mathit{\tau }}_{\mathrm{tot}}\left(\mathit{\lambda }\right)},\end{array}$$

where I_lamp(λ) is the radiance emitted by the lamp and the total optical depth ${\mathit{\tau }}_{tot}\left(\mathit{\lambda }\right)={\mathit{\tau }}_{\mathrm{plume}}\left(\mathit{\lambda }\right)+{\mathit{\tau }}_{\mathrm{BG}}\left(\mathit{\lambda }\right)$ is the combination of the total optical depth of the plume and of the background atmosphere. If an observation is taken when the plume is not in the line of sight of the spectrometer, a background reference spectrum is obtained:

$$\begin{array}{}\text{(2)}& I_{\mathrm{BG}}\left(\mathit{\lambda }\right)=I_{\mathrm{lamp}}\left(\mathit{\lambda }\right)e^{-{\mathit{\tau }}_{\mathrm{BG}}\left(\mathit{\lambda }\right)}.\end{array}$$

By taking two observations with and without the plume in the line of sight, the plume optical depth can be isolated:

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{I\left(\mathit{\lambda }\right)}{I_{\mathrm{BG}}\left(\mathit{\lambda }\right)}}=e^{-{\mathit{\tau }}_{\mathrm{plume}}\left(\mathit{\lambda }\right)}.\end{array}$$

If the plume is considered homogeneous, its optical depth can be expressed as follows:

$$\begin{array}{}\text{(4)}& {\mathit{\tau }}_{\mathrm{plume}}\left(\mathit{\lambda }\right)=\underset{\mathrm{path}}{\int }{k}_{\mathrm{plume}}^{\mathrm{ext}}\left(\mathit{\lambda }\right)\mathrm{d}l=k_{\mathrm{plume}}^{\mathrm{ext}}\left(\mathit{\lambda }\right)L,\end{array}$$

where L is the path length into the plume and $k_{\mathrm{plume}}^{\mathrm{ext}}\left(\mathit{\lambda }\right)$ is the spectral extinction coefficient of the plume, which encompasses both absorption by volcanic gases and absorption and scattering by particles. To avoid radiative interferences with volcanic gaseous effluents other than SO₂, we (1) restricted our analyses to the spectral range 800–1170 cm⁻¹, where only SO₂ (roto-vibrational ν₁ band (e.g. Carboni et al., 2012)) and water vapour have absorption bands, and (2) carefully selected spectral micro-windows within this interval to single out the spectral regions not affected by water vapour lines absorption. The water vapour continuum absorption, in this band, is very small and can be neglected (Shine et al., 2016). In the selected spectral micro-windows, the extinction coefficient of the plume can be expressed as follows:

$$\begin{array}{}\text{(5)}& k_{\mathrm{plume}}^{\mathrm{ext}}\left(\mathit{\lambda }\right)=k_{\mathrm{aer}}^{\mathrm{ext}}\left(\mathit{\lambda }\right)+k_{{\mathrm{SO}}_{\mathrm{2}}}^{\mathrm{abs}}\left(\mathit{\lambda }\right).\end{array}$$

In the previous equation, $k_{\mathrm{aer}}^{\mathrm{ext}}\left(\mathit{\lambda }\right)$ represents the extinction by particulate matter in the plume, i.e. sulfate aerosols, ash or condensed water. Ash emissions were not visually observed during the measurement session. In addition, the spectral signature of ash and pure condensed water, both very different and distinguishable from the SA spectral signature, is not observed in our dataset, so we exclude the presence of both types of particles.

We first selected a background reference spectrum based on the least amount of hydrochloric acid (HCl) measured in the plume, with an independent method, as done by Horrocks et al. (1999). The HCl is abundant in the plume. Then, using Eqs. (3) and (4), we derived the measured $k_{\mathrm{plume}}^{\mathrm{ext}}\left(\mathit{\lambda }\right)$ as follows:

$$\begin{array}{}\text{(6)}& k_{\mathrm{plume},\mathrm{meas}}^{\mathrm{ext}}\left(\mathit{\lambda }\right)={\displaystyle \frac{\mathrm{1}}{L}}\ln {\displaystyle \frac{I_{\mathrm{BG}}\left(\mathit{\lambda }\right)}{I\left(\mathit{\lambda }\right)}}.\end{array}$$

Using an optimal estimation method based on the Levenberg–Marquardt minimisation algorithm (Rodgers, 2000, p. 92–93, and references therein) (more details on the set-up of this method are given in Appendix A), we fitted the measured $k_{\mathrm{plume},\mathrm{meas}}^{\mathrm{ext}}\left(\mathit{\lambda }\right)$ and a modelled $k_{\mathrm{plume},\mathrm{mod}}^{\mathrm{ext}}\left(\mathit{\lambda }\right)$, parameterised as follows:

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& k_{\mathrm{plume},\mathrm{mod}}^{\mathrm{ext}}\left(\mathit{\lambda },\left[{\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}\right]\right)=M_{\mathrm{SA}}^{\mathrm{r}}{k}_{\mathrm{SA}}^{\mathrm{ext}}\left(\mathit{\lambda },\left[{\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}\right],M_{\mathrm{SA}}\right)\\ & +{\mathit{\rho }}_{{\mathrm{SO}}_{\mathrm{2}}}{\mathit{\sigma }}_{{\mathrm{SO}}_{\mathrm{2}}}^{\mathrm{abs}}(\mathit{\lambda },T,p).\end{array}\end{array}$$

In Eq. (7), $k_{\mathrm{SA}}^{\mathrm{ext}}\left(\mathit{\lambda },\left[{\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}\right],M_{\mathrm{SA}}\right)$ is the extinction coefficient of a target sulfate aerosol (SA) layer, calculated using a Mie code (http://eodg.atm.ox.ac.uk/MIE/, last access: 4 October 2019). We fixed the size distribution as a mono-modal log-normal distribution with 0.2 µm mean radius, 30 particles per cubic centimetre number concentration and 1.86 µm width. The SA layers have been taken as a dispersion of spherical droplets of a binary system solution, composed of water and sulfuric acid. The sulfuric acid mixing ratio- ([H₂SO₄]) and temperature-dependent complex refractive indices of these droplets have been taken from Biermann et al. (2000) for the temperature T of the plume, which was assumed to be atmospheric temperature. Each combination of the size distribution parameters and [H₂SO₄] corresponds to a total sulfate aerosol mass concentration. In the spectral region between about 800 and 1200 cm⁻¹, a clear spectral signature of SA has been found by Sellitto and Legras (2016). Sellitto and Legras have also shown that the mass concentration, through a shift of the total absorption signature, and [H₂SO₄], through a modification of the shape of the spectral signature, mostly determine the extinction coefficient of sulfate aerosol layers, except for extreme values of the mean radius (e.g. mean radii greater than 0.4–0.5 µm). In this latter case, the scattering component of the total extinction may become important, leading to a stronger dependence of $k_{\mathrm{plume},\mathrm{mod}}^{\mathrm{ext}}$ on the size distribution. Extreme values of the mean radius are unlikely for these proximal observations. Based on these considerations, in the fitting procedure we adjust the total mass concentration of the sulfate aerosols in the plume by adjusting a mass concentration ratio parameter $M_{\mathrm{SA}}^{\mathrm{r}}$, defined as the ratio between the total mass concentration and the mass concentration associated to the Mie-calculated extinction coefficient $k_{\mathrm{SA}}^{\mathrm{ext}}\left(\mathit{\lambda },\left[{\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}\right],M_{\mathrm{SA}}\right)$. In addition, the optimal estimation is run multiple times, using each time a different value of [H₂SO₄]. An optimal value of [H₂SO₄] is selected based on the deviation of the modelled and measured $k_{\mathrm{plume}}^{\mathrm{ext}}$. As for the SO₂ component, ${\mathit{\sigma }}_{{\mathrm{SO}}_{\mathrm{2}}}^{\mathrm{abs}}(\mathit{\lambda },T,p)$ is the SO₂ absorption cross section calculated at temperature T and pressure p of the plume, using HITRAN (high-resolution transmission molecular absorption database) spectroscopic data (position, shape and intensity of the absorption lines in the selected spectral region). Then, the number concentration ${\mathit{\rho }}_{{\mathrm{SO}}_{\mathrm{2}}}$ is adjusted during the spectral fitting.

From the SA mass concentration ratio $M_{\mathrm{SA}}^{\mathrm{r}}$ and the SO₂ density ${\mathit{\rho }}_{{\mathrm{SO}}_{\mathrm{2}}}$, the total mass concentrations of sulfate aerosols and sulfur dioxide, $M_{\mathrm{SA}}^{\mathrm{T}}$ and $M_{{\mathrm{SO}}_{\mathrm{2}}}^{\mathrm{T}}$, are derived. Finally, the vector of the estimated plume's parameters is tri-dimensional and is composed of the SO₂ and sulfate aerosol mass concentrations, $M_{{\mathrm{SO}}_{\mathrm{2}}}^{\mathrm{T}}$ and $M_{\mathrm{SA}}^{\mathrm{T}}$, and the sulfate aerosol acidity (in terms of their sulfuric acid mixing ratio [H₂SO₄]). The error content of the SO₂ and sulfate aerosol mass concentrations can be derived using the Rodgers framework (Rodgers, 2000). These uncertainties are dominated by error components transferred from the spectral measurement radiometric noise and the information content correlations between these two retrieved parameters, and is typically smaller than 15 % for both parameters.

4 Results and discussion

Figure 1a shows the measured $k_{\mathrm{plume},\mathrm{meas}}^{\mathrm{ext}}$ and modelled (after spectral fitting) extinction coefficient $k_{\mathrm{plume},\mathrm{mod}}^{\mathrm{ext}}$ for the Masaya volcano plume, as well as the spectral residuals, for one individual FTIR observation in our dataset. The raw FTIR observations for this spectrum are in the Supplement. In Fig. 1b, the individual contributions of the SA layer and of the SO₂ to the fitted $k_{\mathrm{plume},\mathrm{mod}}^{\mathrm{ext}}$ are also shown. A clear spectral signature of SA is visible in $k_{\mathrm{plume},\mathrm{meas}}^{\mathrm{ext}}$ with the two peculiar absorption features (two combinations of bend–stretch vibrational modes) of the undissolved H₂SO₄ molecules (870–920 and 1150–1200 cm⁻¹) and the symmetric stretch ν₁ bisulfite ion vibrational mode around 1050 cm⁻¹ (Biermann et al., 2000; Sellitto and Legras, 2016). These components of $k_{\mathrm{plume},\mathrm{mod}}^{\mathrm{ext}}$ are isolated in the SA-only spectrum of Fig. 1b. Comparing the two panel, it is also apparent how an SA-only plume cannot completely explain the spectral shape of the measured spectrum. Adding the SO₂ contribution, the spectrum is more satisfactorily fitted. The residuals are between near-zero, for wavenumbers <1080 cm⁻¹, and about $\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{6}}$ cm⁻¹ (up to 20 %), for longer wavenumbers. The residuals are consistently contained in the ±1σ_residuals interval, except for wavenumbers higher than 1100 cm⁻¹. The higher residuals above 1100 cm⁻¹ can be due to (a) the lack of fine-scale structures of the Biermann et al. (2000) refractive indices for SA (which are provided at a coarse spectral resolution), (b) uncertainties in the SO₂ spectral absorption coefficient transferred from uncertainties in temperature and pressure used in their calculation, and/or (c) spectral drifts in the FTIR observations.

Figure 1 (a) Measured (solid black line) and modelled/fitted (dotted grey line) extinction coefficient of the Masaya volcano plume for the observation of 15 March 1998 at 17:06 UTC. Residuals (modelled minus measured extinction coefficient) are also shown in red. Orange lines show the corresponding standard deviation $\pm \mathrm{1}{\mathit{\sigma }}_{\mathrm{residuals}}=\sqrt{\mathrm{2}}{\mathit{\sigma }}_{k_{\mathrm{plume},\mathrm{meas}}}^{\mathrm{ext}}$ interval (here, excluding specific systematic spectroscopic errors transmitted from the SA laboratory measurements, ${\mathit{\sigma }}_{k_{\mathrm{plume},\mathrm{mod}}^{\mathrm{ext}}}$ has been considered conservatively equal to ${\mathit{\sigma }}_{k_{\mathrm{plume},\mathrm{meas}}^{\mathrm{ext}}}$; details on the calculation of ${\mathit{\sigma }}_{k_{\mathrm{plume},\mathrm{meas}}^{\mathrm{ext}}}$ are given in Appendix A). (b) Individual components of the modelled extinction coefficient due to SA (blue line) and SO₂ (sky blue line).

During the measurement session mentioned in Sect. 2, 38 in-plume spectral measurements are taken. In most cases, i.e. 25 spectra (66 % of the overall observations), coexistent SO₂ and SA are detected. For 5 spectra (13 %) only SO₂ is detected. For the remaining 8 spectra (21 %), no clear SA or SO₂ signals are detected. For these latter cases, the measured plume's extinction coefficient is very small and we suppose that the plume was not or only partially in the line of sight of the FTIR. The mean values of the total mass concentrations of SA and SO₂ and the mixing ratio of SA, calculated by averaging the mentioned 25 individual retrievals obtained during the measurement session, are summarised in Table 1. With an independent retrieval algorithm, using a different spectral range (between 2465 and 2550 cm⁻¹), Horrocks et al. (1999) have found a SO₂ concentration of 98.4±20.0 mg m⁻³ for the same observation set at Masaya. Our retrieved SO₂ average total mass concentration (153.7±85.3 mg m⁻³), even if slightly higher, is consistent with the one obtained by Horrocks et al. (1999). The higher value in our average estimation can be due to the presence of volcanic fluorine compounds that might be present in the plume but are neglected in the present study. The retrieved average sulfate aerosol mass concentration is 0.4±0.2 mg m⁻³. The SO₂∕SA ratio is then about 400. Due to the proximal observations, these aerosols can be considered primary emissions and are still not significantly processed by interaction with the atmosphere. The average H₂SO₄ mixing ratio is 65±18 %. This result reveals highly acidic sulfate aerosols as observed previously with direct sampling (Allen et al., 2002). These acidic aerosols have been shown to be formed by small particles (<2 µm), accounting for up to 80 % of the total mass of aerosol emissions at this volcano. The acidity comes from the rapid conversion of sulfuric acid and hydrofluoric acid to particles, initiated from the oxidation of SO₂ (Mather et al., 2004). Thus, one important source of systematic error in our method is the possible presence of hydrofluoric acid dissolved in the aerosol droplets, that has been neglected in our radiative calculations. By considering the average SO₂ and SA mass concentrations and the average mixing ratio, it results that the gas-to-particle partition of sulfur in the plume is about 400, which is very similar to what was found, in 2002, at Masaya by Allen et al. (2002) (about 450). In any case, the temporal evolution of the three parameters and the SO₂∕SA ratio are very variable and different samples, taken at different times, can be significantly different. Masaya volcanic activity is very stable but the gas emission rates vary and, critically, the aerosol evolution is probably strongly dependent on highly variable atmospheric parameters such as relative humidity and temperature, as well as the plume age at the moment of detection, which depends partly on wind velocity. Notwithstanding this variability, the average sulfur partitioning seems very stable at longer timescales.

Table 1 Mean SO₂ and SA mass concentrations, and mean sulfuric acid mixing ratio in SA, for our retrieval session at Masaya volcano with associated standard deviations (due to their variability during the measurement session). These mean values are obtained by averaging 25 individual retrievals.

We tested the possibility that SA extinction could interfere with SO₂ retrievals by retrieving SO₂ only, with our optimal estimation method, in spectra where the coexistence of the two volcanic plume components is clear, like the one in Fig. 1. In this case, we obtain an overestimation of the SO₂ total mass concentration of nearly 50 % with respect to the case of co-retrieval of SO₂ and SA mass concentrations. Even if this is an extreme case and inversion methods can be developed to partially compensate broadband biases, like the one introduced by aerosols in the sampled air masses, this example shows how the radiative interferences between SO₂ and SA may introduce overestimations in these retrievals. This is particularly important for the observation of volcanic plumes, with a potentially substantial impact in the case of aged plumes. For stratospheric eruptions, rapid formation of secondary SA has been recently observed, e.g. for the Kasatochi volcano (Alaska, USA, August 2008) (Krotkov et al., 2010) and Nabro volcano (Eritrea, June 2011) (Penning de Vries et al., 2014) eruptions, with SO₂ lifetimes as short as a few hours. For volcanic eruptions with injection at lower altitudes (in the troposphere) or in persistent passive degassing regimes, the rapid gas and aqueous-phase oxidation and/or nucleation of SO₂ can lead to radiatively active layers within a few hours after the initial SO₂ emission (Sellitto et al., 2016; Guermazi et al., 2017), with potentially systematic volcanic signature on the regional aerosol optical properties downwind (Sellitto et al., 2017c). In general, the lifetime of SO₂ and the inherent timescales of SA formation and evolution are complex and depend on various factors including the solar irradiation, humidity, temperature, pH and the presence of oxidants (e.g. Eatough et al., 1997). Then, we recommend to either co-retrieve SO₂ and SA or to take explicitly into account these interferences when attempting to retrieve chemically, micro-physically, and radiatively interconnected pollutants. These considerations apply notably to satellite observations where the spectroscopic inversion problem can be even more severe than in the present ground-based OP-FTIR case, due to the strongly ill-posed problem linked to a longer atmospheric path and the inherent radiative transfer. We also highlight that to date the typical band used to quantify SO₂ with ground-based OP-FTIR is that at 2500 cm⁻¹, which is less affected by the aerosol extinction (e.g. Horrocks et al., 1999).

5 Conclusions

We developed a retrieval algorithm to observe volcanic SO₂ and SA emissions using ground-based OP-FTIR spectrometry. To the best of our knowledge, this is the first time that SO₂ and SA are simultaneously characterised using this or other ground-based or satellite-based remote-sensing techniques. The retrieval is based on a non-linear least square fitting algorithm, minimising the difference between modelled and measured spectral extinctions of volcanic plumes. A spectral micro-windows selection was performed in order to avoid the interference with water vapour absorption and to optimise the spectral fitting. The absorption of SO₂ and the extinction of SA have been modelled by means of high-spectral-resolution temperature- and pressure-dependent spectroscopic data and a Mie code driven by state-of-the-art aerosol optical properties, respectively.

We applied this method to proximal ground-based FTIR observations at Masaya volcano. We retrieved SO₂ and SA total mass concentrations and SA mixing ratio (linked to particle acidity). Average concentrations of SO₂, SO₂/SA ratios and acidity of aerosols are consistent with previous observations at Masaya. A gas-to-particle partition of sulfur of about 400 is found and a strong acidity ([H₂SO₄] of about 65 %), which is very consistent with past observations at Masaya, pointing to reasonable long-term stability of Masaya primary sulfate emissions. We also underline the importance of taking into account SA when attempting to retrieve volcanic SO₂ from ground-based or satellite remote-sensing instruments around the 1200 cm⁻¹ SO₂ band. Neglecting the rapid formation (or even primary emission) of SA may lead to significant overestimation of SO₂ (in our case up to about 50  %). This method can be fruitfully applied to the quite large library of previously collected FTIR spectra from volcanoes worldwide to constrain the magmatic and/or atmospheric processes that determine their sulfur emissions and gas-to-particle partitioning. This method can also be exported to thermal infrared (TIR) satellite observations, like those from IASI (Infrared Atmospheric Sounding Interferometer) or others, as long as the full radiative transfer through the longer atmospheric path to the satellite platform is taken into account.

Appendix A: Optimal estimation set-up

An ad hoc optimal estimation retrieval method is used in this work to retrieve $M_{\mathrm{SA}}^{\mathrm{r}}$ and ${\mathit{\rho }}_{{\mathrm{SO}}_{\mathrm{2}}}$ that are used to obtain the final tri-dimensional output vector composed of the SO₂ and sulfate aerosol mass concentration $M_{\mathrm{SO}\mathrm{2}}^{\mathrm{T}}$ and $M_{\mathrm{SA}}^{\mathrm{T}}$ and the sulfate aerosol acidity [H₂SO₄] (see Sect. 3). The method is based on the minimisation of the following cost function J:

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}& J={\left(\mathit{x}-{\mathit{x}}_a\right)}^{\mathrm{T}}{\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}\left(\mathit{x}-{\mathit{x}}_{\mathrm{a}}\right)\\ & +{\left(k_{\mathrm{plume},\mathrm{meas}}^{\mathrm{ext}}-k_{\mathrm{plume},\mathrm{mod}}^{\mathrm{ext}}\left(\mathit{x}\right)\right)}^{\mathrm{T}}{\mathbf{S}}_{\mathit{ϵ}}^{-\mathrm{1}}\\ & \left(k_{\mathrm{plume},\mathrm{meas}}^{\mathrm{ext}}-k_{\mathrm{plume},\mathrm{mod}}^{\mathrm{ext}}\left(\mathit{x}\right)\right).\end{array}\end{array}$$

In the previous equation, $\mathit{x}=\left[M_{\mathrm{SA}}^{\mathrm{r}},{\mathit{\rho }}_{{\mathrm{SO}}_{\mathrm{2}}}\right]$ is the state vector to be retrieved; x_a is the a priori vector for x, with its associated covariance matrix S_a; and S_ϵ is the measurement error covariance matrix (linked to the uncertainty in the measured extinction coefficient $k_{\mathrm{plume},\mathrm{meas}}^{\mathrm{ext}}$). The measured and modelled extinction coefficients $k_{\mathrm{plume},\mathrm{meas}}^{\mathrm{ext}}$ and $k_{\mathrm{plume},\mathrm{mod}}^{\mathrm{ext}}\left(\mathit{x}\right)$ are the same as in Sect. 3, with their implicit wavelength dependence and represented as vector.

The a priori values of $M_{\mathrm{SA}}^{\mathrm{r}}$ and ${\mathit{\rho }}_{{\mathrm{SO}}_{\mathrm{2}}}$, elements of x_a, are both taken as zero. In the a priori covariance matrix S_a, the diagonal elements describe the expected variability in the retrieved parameters, while off-diagonal elements express the possible covariances: in this work only diagonal elements are assigned non-zero values. The measurement noise covariance matrix S_ϵ has been constructed to represent the measurement error on the measured extinction coefficient. Thus, for each diagonal element (measurement noise transferred to the extinction coefficient at a given wavelength), we have considered the following expression, with reference to Eq. (6):

$$\begin{array}{}\text{(A2)}& {\mathit{\sigma }}_{k_{\mathrm{plume},\mathrm{meas}}^{\mathrm{ext}}}={\displaystyle \frac{\mathrm{1}}{L^{\mathrm{2}}}}\left({\displaystyle \frac{{\mathit{\sigma }}_{I_{\mathrm{BG}}}^{\mathrm{2}}}{I_{\mathrm{BG}}^{\mathrm{2}}}}+{\displaystyle \frac{{\mathit{\sigma }}_I^{\mathrm{2}}}{I^{\mathrm{2}}}}\right).\end{array}$$

Based on a conservative choice of a relative error of 1 % for both background and in-plume observations, we obtain a ${\mathit{\sigma }}_{k_{\mathrm{plume},\mathrm{meas}}^{\mathrm{ext}}}$ variance of about $\mathrm{8.0}\times {\mathrm{10}}^{-\mathrm{14}}$ cm⁻². Off-diagonal values of S_ϵ are taken as zero.