The conventional DOAS method applied in routine remote sensing of volcanic SO₂ assumes that all the measured light has been last scattered into the field of view above the plume and has passed in a straight line through the plume. In practice, however, a fraction of the light that reaches the instrument originates from atmospheric scattering on molecules and aerosols below the plume and thus has not passed through the plume. This causes a “dilution” of the slant column retrieved by the DOAS method.

We define the following:

• Sky spectrum: a spectrum of the skylight that has not passed through the plume. This can be measured in a location with plume free zenith view.

Figure 1 shows a theoretical SO₂ transmittance spectrum covering the region commonly used in DOAS evaluations, with the three narrow evaluation intervals used in this paper highlighted (305–310, 310–315, and 315–320 nm). When a certain fraction of the measured light has not passed through the plume, it will dilute the measured spectrum and act to reduce the measured absorption. If we assume that the spectral shape of this diluting light is identical to that of the sky spectrum, then the different absorption bands will be diluted with the same fraction. However, if the absorption is strong enough that the absorption features in the shorter wavelength bands, which have the strongest absorption, do not respond linearly to the dilution, then this dilution of the measured spectrum will change the differential absorptions more in the shorter wavelength bands compared to the longer wavelength bands where the absorption is weaker. Hence, performing a DOAS evaluation of the different bands highlighted above will give different SO₂ column in each band (Fickel and Delgado Granados, 2017).

Also, the SO₂ broadband absorption contributes to this imbalance because although the different bands have equal absolute dilution, expressed as the ratio of the dilution to the sky spectrum, the effective intensity of light coming from the plume is different in the different bands and thus the relative dilution, with respect to the measured spectrum, will differ. Figure 1 shows the absorbance spectrum of 2,000 ppm m SO₂. Here it can be seen that due to the broadband absorption in SO₂ the actual intensity in the plume at the three wavelength bands used is only 27, 45, and 68%, respectively, of the plume intensity without SO₂. Thus, a dilution factor of 0.1 (10% of the normalized sky spectrum contribution to the measured spectrum), will effectively result in a dilution of the differential SO₂ absorption corresponding to about 37, 22, and 15% in the respective wavelength bands. Note that if there is no dilution then this difference in intensity between the evaluation bands caused by the broadband SO₂ absorption is not a problem as it is taken care of by the DOAS evaluation procedure. However, when dilution is present, this difference in intensity contribution from the plume results in a different effective dilution between the bands. In the following we refer to this effect as an “imbalance” in the absorption and make use of this effect to indicate presence of “below-the-plume” scattering and quantitatively correct for its diluting effect. As will be discussed further later in this paper, some other wavelength dependent processes, like plume aerosol extinction and atmospheric ozone concentration, may also contribute to the “imbalance” when dilution is present.

The method presented here relies on the assumption that the fraction of diluting light, i.e., the ratio of the light scattered below the plume, the dilution spectrum, to the spectrum of light scattered above the plume, the sky spectrum, is constant within a narrow spectral interval. By narrow we mean a band not wider than 10 nm.

The method is applied as follows: first, we perform conventional DOAS evaluations on two adjacent or slightly overlapping bands encompassing an interval at least 5 nm and not larger than 10 nm wide (e.g., 305–310 and 310–315 nm) to derive the apparent columns of SO₂ in the two bands, using a measured spectrum outside of the plume as the background. If there is dilution below the plume and the actual SO₂ column is sufficiently high, then the two retrieved columns will differ with the column retrieved at the stronger absorption band being lower than that from the weaker absorption band and both being lower than the true column. Notice that the spectral fit obtained at both bands may still look good (no significant systematics in the fit residual) because the effect of dilution on the spectral shape in such narrow spectral intervals is often not evident. If a difference in columns is found between the two bands, we proceed to iteratively subtract a constant fraction of the sky spectrum from the measured spectrum until a fraction is found where the evaluated columns in the two bands coincide. To obtain a quantifiable parameter for the dilution, the sky spectrum used in the iterative subtraction is normalized to become of equal intensity to the measured spectrum at 350 nm, which is chosen to be close to the SO₂ absorption region but not significantly affected even by strong absorptions of SO₂. Hereby we quantitatively obtain the amount of spectral dilution (the dilution factor k, see above) as well as a corrected value for the column. The process is described in more detail and with illustrations in Figures 2, 3.

In the following we formalize the method. For the sake of clarity, we assume that any instrumental effect that may be present in the plume and sky spectra, such as dark current, offset spectra, or possible changes in the instrumental line function or wavelength calibration, have been corrected.

We can represent the algorithm mathematically in the following way: the measured spectrum, I m λ is expressed as the sum of the plume spectrum, I p λ , and the dilution spectrum, I d λ : I m λ = I p λ + I d λ (1)

It is convenient to relate the measured spectrum to the sky spectrum. Both may be measured under different conditions, for instance at different elevation angles or under unequal background scattering conditions. This results in intensities that may vary even in spectral intervals where absorption by the plume is expected to be negligible. We can start by defining a ‘scaled sky spectrum’ as the product of the sky spectrum by a factor equal to the ratio of the intensities of the measured and sky spectrum at a wavelength with negligible SO₂ absorption. We choose the intensities averaged over 5 nm at 350 nm for this scaling: I S ′ λ = I s λ I m 350   n m I s 350   n m (2)

The scaling is not problematic because it is equivalent to allowing a free scaling between the measured spectrum and the sky spectrum in the standard DOAS formalism.

The dilution spectrum is assumed to be proportional to the scaled sky spectrum I s λ within a narrow wavelength interval: I d λ = k I S ′ λ (3)

With the constant k > 0 defined as the dilution factor.

The plume spectrum must then correspond to the complementary fraction 1 − k of the scaled sky spectrum multiplied by a Beer-Lambert law term with a relative optical depth τ λ between the plume and scaled sky spectra. This relative optical depth not only contains the effect of absorption by SO₂ in the plume but also the effect of all other processes that make both spectra different, namely: difference in O₃ absorption, difference in elastic and inelastic (Ring effect) molecular scattering, or difference in extinction due to aerosols. I p λ = 1 − k I S ′ λ e − τ λ (4)

The meaning of the complementary factor 1 − k is to normalize the transmittance of the measured spectrum with respect to the scaled sky spectrum. We do this by noting that in the limit of zero relative optical depth the measured and scaled sky spectrum must be equal.

Then follows that the measured spectrum is equal to: I m λ = I S ′ λ 1 − k e − τ λ + k (5)

The method we propose iteratively subtracts a fraction of the normalized sky spectrum I s ′ λ from the measured spectrum I m λ , normalizes the new transmittance spectrum to compensate for the subtracted term, and applies the DOAS retrieval on two adjacent wavelength intervals. In each iteration, we find a set of fit parameters that minimize the difference between the logarithms of the modelled and measured spectra. These parameters are the coefficients of a polynomial that represents broadband structures of the relative optical depth (between scaled sky and measured spectra), and the relative column densities of atmospheric gases or pseudo-absorbers (O₃, Ring-effect) with absorption features in the interval of evaluation, as well as an apparent column density of SO₂, S i .

Let us use the notation R i , Δ λ I m λ → S i , Δ λ to represent the retrieval of an apparent column density S i at a spectral interval Δ λ applied to the spectrum I m λ . We apply the retrieval in two adjacent narrow spectral intervals on modified spectra formed by subtracting a fraction of the scaled sky spectrum from the measured spectrum, followed by re-normalization: R i , Δ λ I m λ − x i I s ′ λ 1 − x i → S i , Δ λ (6)

This iteration is applied until we find the constant fraction x i that would make the retrieved columns S i at the two adjacent spectral intervals equal. In other words, we vary x_i until S i , Δ λ 1 = S i , Δ λ 2 . We can see from equations [1–5] that this occurs when x i = k , i.e., when we subtract from the measured radiance the component that was added by scattering of radiation not passing through the plume. In this way, we find not only the true column S , but also the dilution factor k , assumed to be constant along the two wavelength intervals used in the analysis.

Depending on the available amount of UV-light and actual SO₂ column strength, different wavelength intervals can be chosen for the evaluation. The idea is to select two wavelength intervals Δ λ 1 and Δ λ 2 , wide enough to facilitate a good DOAS evaluation (∼5 nm), located close enough to encounter a similar dilution effect from atmospheric scattering, while still having a distinct difference in absorption cross-section. In addition, the absorption cross-section should be large enough so that at least one of the included absorption bands encounters non-linearity due to the Beer-Lambert law, given the actual SO₂ columns. Ideally, the intervals would also be differently affected by broadband SO₂ absorption. These conditions imply that different pairs may be optimal under different conditions; optically weak plumes need to be measured at short wavelengths to obtain enough non-linearity effects in the absorption, while optically thicker plumes, may need to be measured at longer wavelengths and thus weaker absorption as the stronger absorption bands may be saturated.

For the method to be able to determine a corrected column and dilution factor it is necessary that the difference between the columns derived using the two different bands can be clearly distinguished so that a well-defined cross-over point can be established. Thus, the true column and dilution must be high enough that the Beer-Lambert law non-linearity or SO₂ broad-band absorption cause a difference between the columns in the bands used that overcomes the noise level. This is illustrated in Figure 4, which shows how the columns derived in the three evaluation bands: 305–310, 310–315, and 315–320 nm are affected by dilution, for SO₂ columns of 500, 1,000 and 2000 ppm m, respectively. These plots are based on theoretical data. In a real case a difference between the columns of 5–10 ppm m is needed to overcome the noise level.

For SO₂ columns below 500 ppm m, the columns derived from the two lower bands, 305–310/310–315, respond almost linearly to increased dilution and the method will not converge and thus will fail to quantify the dilution. However, the dilution is still there and affects the column data linearly (15% dilution gives approximately 15% lower SO₂ column). For SO₂ columns higher than about 600 ppm m the difference between the two lower bands, at 10% dilution, is more than 10 ppm m and thus this pair can be used for evaluation of the dilution. In the same way, for SO₂ columns higher than about 800 ppm m also the higher pair 310–315/315–320 gives an imbalance response enough for evaluation at 10% dilution. At the other end, for SO₂ columns higher than 5,000 ppm m, saturation effects due to the limited spectral resolution start to cause problems for the lower pair. The longer-wavelengths pair is for the same reason limited to about 8,000 ppm m. Thus, for low to moderate SO₂ columns (600–2,000 ppm m) the lower pair is preferred due to its higher sensitivity, provided that the UV light level is high enough at 305 nm to overcome the noise level. For medium to high SO₂ columns (1,000–5,000 ppm m) the higher pair may instead be used, especially when the UV level at 305–310 nm is too low for accurate measurements. Of course, it is not possible to know the actual column before the correction, therefore evaluation at three bands may be necessary to explore the range of variation before choosing the pair of wavelength intervals that would maximize the quality of the retrieval based on the considerations given above.

In the examples and discussion above we assumed a dilution of 10%. With less dilution even higher true SO₂ columns are needed to give a reliable detection. This is also illustrated in Figure 4. The band with weakest absorption (315–320 nm) drops almost linearly with increasing dilution, while the band with the strongest absorption drops off faster than linearly initially and then drops more linearly, almost in parallel to the weak band. This is because with increasing dilution, the apparent column gets smaller and smaller and eventually the stronger absorption bands also become linear. By following the two curves belonging to a specific true SO₂ column from high to low dilution, one can simulate how the dilution correction method works. Notably, in the case of strong dilution (for example 50%) and strong true SO₂ column (for example, 2,000 ppm m), the apparent columns in the two stronger absorption bands initially diverge when dilution is decreased and then, as dilution continues to decrease, the lines of retrieved columns at the two bands start to converge. Thus, the iterative correction algorithm must allow for both convergence and divergence until the retrieved values are close enough and the column obtained at the weaker band is not lower than that obtained at the stronger band.

It is also interesting to note that, although the two processes of non-linearity and broad-band absorption both act to enhance the imbalance described above, the two processes in a way counterbalance each other. Figure 5 shows how the apparent SO₂ column, evaluated in the evaluation window 309–314 nm, varies with true SO₂ column for different dilution factors.

For a given dilution constant, the effect of dilution due to broad-band SO₂ absorption becomes stronger and stronger with increasing SO₂ column density at the same time as the monotonous increase in apparent column due to non-linearity gets smaller and smaller. This results in a turnover making the apparent column decrease with increasing true SO₂ column. As a consequence, for a given evaluation band and dilution factor there is a maximum column that can be obtained in a DOAS evaluation. This effect has been pointed out by others (Mori et al., 2006; Varnam et al., 2020) and constitutes a limitation in standard DOAS evaluations as well as in some existing schemes for dilution correction. However, this is not a limitation for our method, as we just keep iterating until we reach the convergence point.

Besides the effect of dilution, conditions where light with dramatically different SO₂ column are superimposed in a single spectrum may also cause such an imbalance. This can occur when spectra containing different columns are co-added into a single spectrum, or if the field-of-view of the receiving optics cover a wider angle than represented by typical variations in total column in the plume, e.g., when a single spectrum contains light from both strong and weak columns. A typical example of this situation is when measurements are made at the beginning and end of a traverse where dramatic changes in column occur over short traverse distances. These effects can, to some extent, be avoided by using a narrow field-of-view and high sampling rate, and by excluding data from the traverse edges in the scattering correction dilution estimate.

Other processes resulting in similar effects are straylight in the spectrometer and scattering within the plume itself. The amount of straylight in the spectrometer depends on type of spectrometer as well as its condition and can to some extent be reduced by filtering out light with wavelengths outside the studied wavelength region. Moreover, the impact that straylight has on the imbalance can be checked by performing a measurement with a gas cell using a clean air background, with high enough gas column in the cell to reveal imbalance in the extinction law, with a nominal column value of 1,000 ppm m or more. If an imbalance in the absorption is detected in this cell measurement, straylight is likely to be present.

Eliminating effects of scattering within the plume is more difficult. One approach to study the impact of scattering within the plume on the spectral imbalance used in our method is modelling. The results from a study using the radiative transfer model McArtim (Deutschmann et al., 2011) is shown in the chapter 5.1 below. This modelling study indicates that the method presented here may give somewhat reduced column densities, especially for the case of a non-condensing plume above an opaque atmosphere. Plumes with high aerosol optical depths tend to offset this effect producing higher columns, probably due to multiple scattering.

A crucial condition in the method is that the clean air spectrum used in the iterative subtraction, and the true dilution spectrum caused by scattering below the plume, have similar shape over the two wavelength intervals used in the analysis. If this condition is not fulfilled, then the “crossover point” between the bands in the iterative process will be wrong. This problem is minimized by using two narrow wavelength bands located close to each other and to perform the sky and plume measurements close in time. Special care should be taken when cloud conditions differ significantly between sky and plume measurements. One way to check this is to compare the slope of the spectra in the actual wavelength interval, after dividing sky spectra taken before and after the plume traverse. Ozone has strong broadband absorption in the spectral regions suggested here and a difference in ozone between the sky spectrum and the plume spectrum may have an impact on the imbalance between the two evaluation bands. However, the impact of ozone is dominated by its stratospheric slant column density (O₃SCD), which is determined by the solar zenith angle (SZA, Eq. 6). Thus, if the sky spectrum is taken close in time to the plume spectrum, or if measurements are made close to noon, then the impact of ozone is negligible. If these conditions are not fulfilled, then the effect of ozone may still be eliminated by evaluating the relative ozone slant column density before and after passing the plume in the traverse. Using GPS data, time, and location, one can then derive the SZA related to each of these measurements and calculate the ozone vertical column density (O₃VCD) for the time of the traverse. Knowing the O₃VCD one can then use the GPS data of each spectrum to derive the O₃SCD of each spectrum and use the known absolute ozone absorption cross-sections to eliminate the effect of the difference in O₃SCD between the sky spectrum and the plume spectrum (including the dilution spectrum).

If we denote:

O 3 S C D s k y ozone slant column density in sky spectrum

O 3 S C D p l u m e ozone slant column density in plume spectrum

O 3 V C D ozone vertical column density (common to sky and plume)

S Z A solar zenith angle (different for sky and plume) O 3 V C D = O 3 S C D p l u m e − O 3 S C D s k y 1 cos S Z A p l u m e − 1 cos S Z A s k y (6a)

The method described above has been implemented in software, ScatteringCorrection, written in C++ and the initialization parameters are further described in the Supplementary Annex to this article. The software reads measured spectra saved in the extended standard (*.std) format, e.g., used by the MobileDOAS software (Johansson, 2009) and, for each measured spectrum, performs DOAS evaluations in two aforementioned narrow wavelength intervals and in one wider wavelength interval. If the SO₂ column in the wide wavelength interval exceeds a user-defined limit (approx. 500 ppm m), determined from when the non-linearity of the Beer-Lambert law begins, and there is a difference in SO₂ column between the two short wavelength intervals, then the software will attempt to perform a dilution correction of the measurement as described above.

The software performs a search for the optimal dilution factor k by attempting to minimize the difference between the evaluated SO₂ columns in the bands of longer and shorter wavelengths λ₁ and λ₂. S i , Δ λ 1 < S i , Δ λ 2 implies the applied dilution factor k is too small, while S i , Δ λ 1 > S i , Δ λ 2 implies k is too large. This fact is used to determine the direction of the search in each step of the iteration in which k is varied. Notice, however, that if the search has stepped too far beyond the optimum dilution factor k, the above relationship may fail. This fact must be accounted for in the search, most easily by also verifying that the DOAS fit is valid, e.g., by requiring a maximum permissible standard deviation in the fit residual. In our implementation, the search begins by attempting to determine the upper bound for k by stepping up k in relatively large steps until S I , Δ λ 1 > S i , Δ λ 2 . During this stepping up of k, the software also determines the largest lower bound on the dilution factor, as the largest k_i, which yielded S i , Δ λ 1 < S i , Δ λ 2 . The second phase of the search is to determine the optimum dilution factor k between the upper and lower bounds. In the case of valid DOAS fits for both the lower and upper bounds, this is done using a linear interpolation between the two as this method has been found to be the fastest. In the case of an invalid fit on either bound, which typically occurs when the upper bound has stepped too far, a binary search is applied instead.

To improve our understanding on how different parameters like plume height, SO₂ concentration, and aerosol burden in and below volcanic plumes influence our light dilution correction, we applied the method to spectra simulated with the radiative transfer model McArtim (Deutschmann et al., 2011). The model can be used to simulate UV/vis spectra expected in a wide range of user-defined atmospheric and plume conditions, and thus allows us to test the ability to derive accurate SO₂ SCDs using our scattering correction technique. For the purposes of this study, the spectra produced by the model were adapted to standard DOAS std format and then analyzed using the developed ScatteringCorrection program (see details in Supplementary Annex).

We study three scenarios: 1) a clean Rayleigh atmosphere and a plume with zero aerosol optical depth (clean atmosphere/clean plume); 2) the same plume but in an opaque atmosphere where a 1 km thick layer with aerosol optical depth of 0.2 is placed at ground level (opaque atmosphere/clean plume); and 3) an opaque plume with aerosol optical depth of 1, immersed in a clean atmosphere (clean atmosphere/opaque plume). For the three scenarios the plume has a vertical column density of 2,000 ppm m and has a cylindrical shape with 500 m diameter and horizontal orientation. The plume altitude is varied between 0 and 3,000 m above ground level (plume-bottom), which is assumed to be at 50 m above sea level. The Sun shines with a relative azimuth of 90° with respect to the plume axis and a solar zenith angle of 30°. The field of view of the instrument is set to 8 mrad. The instrument is located below the plume axis at 50 m altitude and pointed to the zenith. Other parameters of the simulation are given in the Supplementary Appendix.

Dividing each spectrum of the three scenarios by the respective clean air “sky” spectrum produces transmittance spectra, as shown in Figure 6.

Some correction algorithms have been proposed on the assumption that the diluting light follows the Koschmieder formula using either a Rayleigh scattering coefficient, or another coefficient suggested by the measurement conditions (Koschmieder, 1938; Mori et al., 2006; Vogel et al., 2011; Campion et al., 2015; Varnam et al., 2020). One way to test the assumptions of our model, that assumes a wavelength-independent dilution term, and the predictions of a wavelength-dependent term, is to test predictions of both on synthetic spectra. Our proposed model assumes that the dilution spectrum is a constant fraction of the sky spectrum, whereas in other methods the dilution spectrum is equal to the sky spectrum scaled by a wavelength dependent scattering efficiency. Using the spectra simulated with McArtim in which the true SO₂ column density is known we derived the wavelength-dependency of the dilution factor k (Figure 7).

Our results indicate that the trend in the dilution term is nearly independent of wavelength, a circumstance that is in good agreement with our proposed model. There are some structures that are correlated with the absorption of SO₂ which may be due to incomplete compensation of the term involving SO₂ in our model. This may be due to radiative transfer effects that produced an effective column that is not equal to the prescribed vertical column density used to model the different spectra. In the interpretation given by Kern et al. (2010), this could be caused by different effective pathlengths due to different absorption strength of SO₂ for different wavelengths.

Applying the ScatteringCorrection software using the lower bands 305–310/310–315 nm, we obtain the following results for retrieved column and dilution factor (Figure 8).

The modelling results indicate that our proposed method can deal relatively well with spectra taken under conditions of clear and opaque atmospheres and clear and opaque plumes. In all scenarios, we notice the increase of dilution factor with distance to the plume and with increase in scattering conditions below the plume. The retrieved SO₂ column densities without correction show a systematic decrease with increasing dilution, but aerosols in the plume seem to reduce the amount of dilution, possibly due to enhanced forward scattering towards the instrument. The retrieved column values after correction lie within 10% accuracy of the nominal values, being worst for the case of an opaque atmosphere and slightly better for the case of an opaque plume. We think this reflects the effect of two processes: one is that the percentage of radiation that reaches the instrument passing through the plume may be larger for an opaque plume due to higher forward scattering, and the other is that the aerosols may increase the mean path of radiation in the plume due to multiple scattering. The modelling results also show that the assumption of constant dilution within a narrow wavelength interval is reasonable.

Next, we tested our scattering correction on real-world DOAS traverse measurements. Three different cases are shown, representing different column densities and scattering conditions: Fagradalsfjall located on Reykanaes peninsula on south Island (low columns 800 ppm m, clean air), Manam located off the northeast coast of mainland Papua New Guinea (very high columns 5,000 ppm m, high scattering), and Holuhraun located in the inland of north east Island (moderate columns 1,000–2,000 ppm m, low scattering).