Improved optical flow velocity analysis in SO 2 camera images of volcanic plumes – implications for emission-rate retrievals investigated at Mt Etna , Italy and Guallatiri , Chile

Accurate gas velocity measurements in emission plumes are highly desirable for various atmospheric remote sensing applications. The imaging technique of UV SO2 cameras is commonly used to monitor SO2 emissions from volcanoes and anthropogenic sources (e.g. power plants, ships). The camera systems capture the emission plumes at high spatial and temporal resolution. This allows the gas velocities in the plume to be retrieved directly from the images. The latter can be measured at a pixel level using optical flow (OF) algorithms. This is particularly advantageous under turbulent plume conditions. However, OF algorithms intrinsically rely on contrast in the images and often fail to detect motion in low-contrast image areas. We present a new method to identify ill-constrained OF motion vectors and replace them using the local average velocity vector. The latter is derived based on histograms of the retrieved OF motion fields. The new method is applied to two example data sets recorded at Mt Etna (Italy) and Guallatiri (Chile). We show that in many cases, the uncorrected OF yields significantly underestimated SO2 emission rates. We further show that our proposed correction can account for this and that it significantly improves the reliability of optical-flow-based gas velocity retrievals. In the case of Mt Etna, the SO2 emissions of the northeastern crater are investigated. The corrected SO2 emission rates range between 4.8 and 10.7 kgs−1 (average of 7.1 ± 1.3kgs−1) and are in good agreement with previously reported values. For the Guallatiri data, the emissions of the central crater and a fumarolic field are investigated. The retrieved SO2 emission rates are between 0.5 and 2.9 kgs−1 (average of 1.3 ± 0.5kgs−1) and provide the first report of SO2 emissions from this remotely located and inaccessible volcano.


Introduction
Studying and monitoring gas emissions is highly desirable since the emitted gases can have substantial environmental impacts.This includes both natural and anthropogenic sources such as volcanoes, industrial areas, power plants, urban emissions or wildfires.The measurements can help to better assess regional and global impacts of the emissions, for instance, related to air-quality standards and pollution monitoring or climate impacts (e.g.Schwartz, 1994, Robock, 2000, IPCC, 2013).
Sulfur dioxide (SO 2 ), in particular, is a toxic gas emitted both by anthropogenic and natural sources (e.g.power plants, ships, volcanoes).The pollutant has various impacts, both of socio-environmental and economic nature (e.g.human health, agriculture) and on the climate (e.g.being a precursor of stratospheric sulfur aerosols, Wigley, 1989).Furthermore, SO 2 is an important monitoring parameter related to volcanic risk assessment (e.g.Fischer et al., 1994, Caltabiano et al., 1994).
Passive remote sensing techniques are commonly used to monitor gas emissions from localised emitters (or point sources).The instruments are based on the principle of light absorption and typically measure path-integrated concentrations (column densities, CDs) of the gases.Instrumentation can be ground, airborne and satellite based and can cover wavelengths ranging from the near ultraviolet (UV) up to thermal long-wave infrared (LWIR), either using solar or thermal radiation as a light source.Note that the term "point source" is not clearly defined and may, in some cases, refer to scales of several kilometres (e.g. a whole city in the case of space-based observations), and in other cases, to only a few metres (e.g. a power-plant chimney for ground-based near-source measurements).
Gas emission rates (or fluxes) of the sources are typically retrieved along a plume transect by integrating the product of the measured CDs with the local gas velocities in the plume.The latter may be estimated using meteorological weather data (e.g.Frins et al., 2011) or using correlation techniques (e.g.Williams-Jones et al., 2006, Johansson et al., 2009) if the measurements are performed at a moderate sampling rate (e.g.spectroscopic instrumentation such as COSPEC or scanning DOAS instruments, e.g.Moffat andMillan, 1971, Platt andStutz, 2008) and at sufficient source distance.
A more recent measurement technique is based on camera systems which are equipped with wavelength selective filters (e.g.McElhoe and Conner, 1986, Mori and Burton, 2006, Prata and Bernardo, 2014, Kuhn et al., 2014, Dekemper et al., 2016).The imaging devices can be used to create instantaneous CD maps of the measured species (e.g.SO 2 , NO 2 ) at high spatial resolution and at sampling rates potentially down into the sub-Hz regime (depending on the optical set-up and lighting conditions).This allows us to study high-frequency variations in the emission signals or to investigate individual sources separately (e.g.Tamburello et al., 2013, D'Aleo et al., 2016).As a result, the cameras are often pointed at the vicinity of the source, where the plumes can show turbulent behaviour, mostly as a result of aerodynamic effects and buoyancy.The resulting velocity fields often deviate significantly from the meteorological background wind field.Luckily, the high resolution of the imaging systems allows us to account for these spatial and temporal fluctuations by directly measuring the projected 2-D velocity fields using optical flow (OF) algorithms (e.g.Krueger et al., 2013, Bjornsson et al., 2013, Peters et al., 2015, Lopez et al., 2015, Stebel et al., 2015, Kern et al., 2015).
OF algorithms can detect motion at the pixel level by tracking distinct image features in consecutive frames.In the following, the basic principles of the OF computation are briefly introduced as well as different optimisation strategies (see e.g.Jähne, 1997, Fleet and Weiss, 2006, Fortun et al., 2015 for a comprehensive introduction into the topic).OF algorithms are based on the assumption that a certain image quantity, such as the brightness I or the local phase φ, is conserved between consecutive frames.Then, a continuity equation of the form (1) can be used to describe the apparent motion of brightness (or phase) patterns between two frames.Here, f = [u, v] T denotes the flow vector in the detector coordinate system i, j .g is the conserved quantity (e.g.I, φ), ∇ ij = [∂ i , ∂ j ] T and ∂ t denote the spatial and temporal differentiation operators.Equation ( 1) is typically referred to as the optical flow constraint (OFC) equation and can be solved numerically per image pixel, for example using a least-squares or a total leastsquares optimisation scheme.The OFC states an ill-posed problem, as it seeks to find the two velocity components u and v from a single constraint (i.e.I or φ; cf.Eq. 1).This is commonly referred to as the aperture problem and is typically accounted for by introducing further constraints that impose spatial coherency to the flow field.These can be subdivided into local and global constraints or a combination of both (e.g.Bruhn et al., 2005).Local methods (e.g.Lucas and Kanade, 1981) apply the coherency constraint only within a certain neighbourhood around each pixel (the size of this aperture can usually be set by the user).Thus, for pixel positions that do not contain at least one trackable feature within the neighbourhood specified by the aperture size, the algorithm will fail to detect motion.We shall see below that this can be a fundamental problem for the emission-rate analysis using plume imagery, in the case that extended homogeneous plume regions coincide with a retrieval transect .The problem is less pronounced for OF algorithms using global constraints (e.g. the algorithm by Horn and Schunck, 1981 which is used in Kern et al., 2015), which can propagate reliable motion vectors over larger image areas.However, note that, depending on the optimisation strategy, global regularisers are often more sensitive to noise (e.g.Barron et al., 1994) and are typically computationally more demanding (e.g.Fleet and Weiss, 2006).
Most of the modern OF algorithms include a multi-scale analysis where the flow field is retrieved from coarse to fine features using image pyramids combined with suitable warping techniques (e.g.Anandan, 1989).This can significantly increase the robustness of the results and is of particular relevance in the case of large displacements (i.e.several image pixels, e.g.Beauchemin andBarron, 1995, Fleet andWeiss, 2006).
Optical flow intercomparison benchmarks (e.g.Baker et al., 2011, Menze andGeiger, 2015) can provide useful information with which to assess the performance (e.g.accuracy) and applicability (e.g.computational demands, availability of source code) of different OF algorithms.Particularly important for the emission-rate analysis is the computational efficiency as well as the performance within homogeneous image regions.As discussed above, the latter may be optimised via the incorporated coherency constraints (e.g. by increasing the local averaging neighbourhood around each pixel) or by performing a multi-scale analysis.However, this can lead to a significant increase in the required computation times (e.g.Fleet and Weiss, 2006) and may therefore be inapplicable, especially for near-real time analyses.
Given these challenges, in many cases the choice of a suitable OF algorithm will be a trade-off between computational efficiency and the performance within homogeneous image regions.In order to rule out potential failures in the OF retrieval, it is therefore highly desirable to assess the OF performance before calculating the emission rates.In this paper, we propose a new method, which analyses an OF displacement vector field (DVF) in order to identify and correct for potentially unphysical OF motion estimates.The correction is performed in a localised manner, within a specific region of interest (ROI) in the images (e.g. in proximity to a plume transect ).It measures the local average velocity vector (LAVV) within the ROI, based on distinct peaks in histograms computed from the local DVF.The strengths of the method are (1) that it is independent of the choice of the OF algorithm and (2) that the additional computational demands are small compared to the OF computation time.The new method is introduced using the Farnebäck optical flow algorithm (Farnebäck, 2003) which showed promising results in Peters et al. (2015) and which is freely available in the OpenCV library (e.g.Bradski, 2000).We use two different volcanic data sets recorded at Mt Etna, Italy and Guallatiri, Chile to show that our method can successfully detect and correct for unphysical OF motion estimates during the emission-rate analysis.
The paper is organised as follows: Sect. 2 starts with a short introduction into the technique of UV SO 2 cameras and the required data analysis.Section 2.2 provides information about the two data sets (i.e.technical set-up, measurement locations), followed by details regarding the image analysis of both data sets (Sect.2.3).The proposed correction for optical-flow-based velocity retrievals is introduced in Sect.2.4.In Sect. 3 the retrieved SO 2 emission rates for the Etna and Guallatiri data sets are presented and compared to results based on (1) the uncorrected OF DVF and (2) assuming a constant global plume velocity using the crosscorrelation lag of integrated plume intersections (e.g.McGonigle et al., 2005).A summary and discussion is given in Sect.4, followed by our conclusions.

UV SO 2 cameras
UV SO 2 cameras measure plume optical densities (ODs) in two wavelength windows of about 10 nm width using dichroic filters.The two filters are typically centred around 310 nm (SO 2 "on-band" filter, i.e. sensitive to SO 2 absorption) and, at nearby wavelengths, around 330 nm (SO 2 "offband" filter).The latter is used for a first-order correction of aerosol scattering in the plume (e.g.Kern et al., 2010).An apparent absorbance (AA) of SO 2 can then be calculated based on the ODs measured in both channels: Here, I , I 0 denote the measured plume and corresponding background intensities.Note that all quantities in Eq. ( 2) are a function of the detector pixel position i, j (e.g.τ AA → τ AA (i, j )).The calibration of the measured AA values (i.e.conversion into SO 2 column densities S SO 2 (i, j )) can be performed using SO 2 calibration cells or using data from a DOAS spectrometer viewing the plume (Lübcke et al., 2013) or a combination of both.The SO 2 emission rates are typically calculated along a suitable plume cross section (PCS) in the SO 2 -CD images S SO 2 (i, j ) (e.g. a straight line) by performing a discrete integration of the form: where m denotes interpolated image coordinates (i, j ) along , f is the camera focal length, d pl is the distance between the camera and the plume and s(m) is the integration step length (for details see Gliß et al., 2017).The effective velocity is measured relative to the normal n of (i.e.constant in the case of straight retrieval lines) using the corresponding velocity vector v(m).The velocities, if retrieved from the images, represent averages along the line of sight (LoS) of each pixel (see e.g.Krueger et al., 2013 for a derivation).Since the velocity components in the LoS direction cannot be measured from the images, the measured velocities are approximately underestimated by a factor of cos(α) (α being the angle between plume direction and image plane).However, to first order (and at small angles α), this cancels out since the length of the LoS inside the plume (and thus, the measured SO 2 -CDs) increases by approximately the same cos(α) factor (Mori and Burton, 2006).

Example data
The proposed method used to correct for unphysical OF velocity vectors is applied to two volcanic data sets recorded at Mt Etna (Italy) and Guallatiri volcano (Chile).Both data sets were recorded using a filter-wheel-based UV SO 2 camera including a DOAS spectrometer.Details about the technical set-up for both data sets are summarised in Table 1.

Etna data
Mt Etna is a stratovolcano situated in the eastern part of the island of Sicily, Italy.We present a short UV camera   1 for a technical set-up of the instruments used for the observations).The data were recorded during a field campaign which took place about 2.5 months prior to a major eruptive event (i.e. in early December 2015, e.g.Smithsonian-Institution, 2013a).The volcano showed quiescent degassing behaviour during all days of the campaign.The measurements were performed from the rooftop of a building located in the town of Milo, about 10.3 km from the source.An overview map is shown in Fig. 1.

Plume conditions
During the 15 min of data, the meteorological conditions were stable, showing a slightly convective plume of the Etna north-eastern crater (NEC) advected downwind (into the left image half; cf.Fig. 1).The emissions of the other craters are more diffuse and could not be fully captured since they were partly covered by the volcanic flank.Therefore, we kept the focus on the NEC emissions, which were investigated along two example PCS lines located at two different positions downwind of the source (orange and blue lines in Fig. 1).
A video of the Etna emissions is shown in the Supplement video no. 1.

Guallatiri data
Guallatiri (18 • 25 00 S, 69 • 5 30 W, 6.071 m a.s.l.) is a stratovolcano located in the Altiplano, northern Chile.The last confirmed eruptive events date back to 1960 (Smithsonian-Institution, 2013b).Due to its remote location little is known about the volcano.The presented data are part of a short field campaign between 20 and 22 November 2014.During the 3 days, the volcano showed quiescent degassing behaviour from the central crater and from a fumarolic field on the SW flank of the volcano.Due to frequent cloud abundances, only a small fraction of the acquired data was suited for the investigation of the SO 2 emissions.A cloud-free time window between 14:48 and 14:59 UTC on 22 November 2014 was chosen (see Ta- ble 1 for details about the instrumental set-up).An overview map is shown in Fig. 2. The measurements were performed at a distance of 13.3 km away from the source.

Plume conditions
Compared to Etna, the Guallatiri emissions showed rather turbulent behaviour with strong variations in the local velocities.The central crater plume, in particular, changed its overall direction significantly over time, which can be seen in the Supplement video no. 2. Emission rates were retrieved along two (connected) PCS lines in the young plume shown in Fig. 2. The lines were chosen such that the emissions from the central crater and the fumarolic field could be investigated separately.

Data analysis
The image analysis was performed using the Python software Pyplis (Gliß et al., 2017).In a first step, all images were corrected for electronic offset and dark current followed by a first-order correction for the signal dilution effect.The latter was applied based on Campion et al. (2015) using suitable volcanic terrain features in the images to retrieve an estimate of the atmospheric scattering extinction coefficients in the viewing direction of the camera.The extinction coefficients were used to correct the measured radiances of plume image pixels for the scattering contribution.The latter were identified using an appropriate τ threshold applied to on-band OD images.
The sky background intensities (required to for the retrieval of AA images, Eq. 2) were determined using on/off sky reference images (SRI) recorded close in time to the plume image data.The background retrieval was done using the background modelling methods 6 (Etna) and 4 (Guallatiri) of the used analysis software Pyplis (Gliß et al., 2017;cf. Table 2 therein).Variations in the sky background intensities and curvature between the plume images and SRI were corrected both in the horizontal and vertical directions using suitable gas (and cloud)-free sky reference areas in the plume images.All AA images were corrected for cross detector variations in the SO 2 sensitivity using a correction mask calculated from cell calibration data as outlined by Lübcke et al. (2013).The AA images were calibrated using plume SO 2 -CDs retrieved from a co-located DOAS instrument (cf.Table 1; see Sect.B1 for details regarding the DOAS retrieval).The position and extent of the DOAS-FOV (field of view) within the camera images are shown in Figs. 1 (Etna) and 2 (Guallatiri) and were identified using the Pearson correlation method described in Gliß et al. (2017).
The gas velocities in the plume were retrieved both using the Farnebäck OF algorithm and the cross-correlation method outlined in McGonigle et al. (2005).Non-physical OF motion vectors along the emission-rate retrieval lines were identified and corrected for using the proposed OF histogram method, which is described in Sect.2.4.Note that for the analysis all images were downscaled by a factor of 2 (using a Gaussian pyramid approach).

Etna
The required plume distances for the emission-rate retrieval were derived from the camera location and viewing direction and assuming a meteorological wind direction of (0 ± 20) • (north wind; cf.Fig. 1).The latter was estimated based on visual observation.The camera viewing direction was retrieved using the position of the south-eastern (SE) crater in the images.The signal dilution correction was performed using atmospheric scattering extinction coefficients retrieved 20 min prior to the presented observations (i.e. from one on and one off-band image recorded at 06:45 UTC; cf.Fig. 10 in Gliß et al., 2017).During this time the camera was pointed at a lower elevation angle and the images contained more suitable terrain features for the correction.Extinction coefficients of on = 0.0743 km −1 and off = 0.0654 km −1 were retrieved and used to correct plume image pixels.The latter were identified from on-band OD images using a threshold of τ on = 0.05.The dilution-corrected AA images were calibrated using the DOAS calibration curve shown in Appendix B2.The linear calibration polynomial was retrieved prior to the analysis using camera AA values that were not corrected for the signal-dilution effect and the corresponding SO 2 -CDs measured with the DOAS spectrometer (for details see Appendix B2).

Guallatiri
The plume distances were retrieved per pixel column assuming a meteorological wind direction of (320 ± 15) • .The latter was estimated based on visual observation combined with a MODIS image (see Supplement) recorded at 15:05 UTC, in which the plume was identified.The viewing direction of the camera was retrieved based on the geographical location of the summit area in the images.
The dilution correction was performed using scattering extinction coefficients of on = 0.0855 ± 0.0012 km −1 and off = 0.0710 ± 0.0008 km −1 .The latter were retrieved between 14:48 and 14:59 UTC using images from a second UV camera, which was equipped with a f = 25 mm lens (i.e. a wider FOV) and hence contained more suitable topographic features for the retrieval.Plume pixels for the dilution correction were identified from on-band OD images using a threshold of τ on = 0.02.An example dilution-corrected SO 2 -CD image is shown in Fig. 2. The DOAS calibration curve is shown in Appendix B2.

Radiative transfer effects
Both the Etna and Guallatiri data were recorded at long distances (> 10km).Consequently, the applied dilution correction accompanies relatively large uncertainties of statistical nature, which we estimate to ±50 %, based on Campion et al. (2015).Furthermore, in-plume radiative transfer (e.g.multiple scattering due to aerosols, SO 2 saturation; see e.g.Kern et al., 2013) may have affected the results to a certain degree.However, both plumes showed only little to no condensation.We therefore assess the impact of aerosol multiple scattering to be negligible.In the case of Etna, SO 2 saturation around 310 nm may induce a small systematic underestimation in the SO 2 emission rates.This is due to the compara-tively large observed SO 2 -CDs of up to 5 × 10 18 cm −2 .The impact of SO 2 saturation is, however likely compensated to a certain degree, since the DOAS SO 2 -CDs (used to calibrate the camera) were retrieved at less affected wavelengths between λ 0 ≈ (315-326) nm (cf.Appendix B1).The same fit interval is used in Gliß et al. (2015), who performed MAX-DOAS measurements of the Etna plume under comparable conditions.They account for SO 2 saturation by using the weak SO 2 bands between λ 1 ≈ (350-373) nm (see also Bobrowski et al., 2010) and find relative deviations of about 10 % between the two wavelength ranges and for SO 2 -CDs exceeding 5 × 10 18 cm −2 (i.e.λ 0 λ 1 ≈ 0.9 cf.Fig. A3 in Gliß et al., 2015).We therefore estimate the impact of SO 2 saturation to be below 20 % for our data.

Optical flow histogram analysis
We developed a method to improve OF-based gas velocity retrievals needed for the analysis of SO 2 emission rates (Eq.3) using UV camera systems.The OF analysis of an image pair yields dense displacement vector fields (DVFs) of the observed gas plumes.In some areas of the image, the DVF represents the actual physical motion of gas in the plume, while other image areas may contain unphysical motion vectors (e.g. in low-contrast plume regions; cf.Sect.1).The proposed method aims to identify all successfully constrained motion vectors and, from these, derives an estimate of the average (or predominant) velocity vector in the plume.The latter is then used to replace unphysical motion vectors in the DVF.We recommend performing the analysis in a localised manner, within a specific region of interest (ROI) since the velocity fields can show large fluctuations over the entire image (e.g.change in direction or magnitude).Based on these histogram peaks, the proposed method derives the local predominant displacement vector (PDV) |f |.A detailed mathematical description of the analysis is provided in Appendix A. In the following, the most important steps of the analysis are described.
The retrieval of the PDV starts with a peak analysis of M ϕ and investigates whether a distinct and unambiguous peak can be identified in the histogram.If this is the case, the expectation value for the local movement direction ϕ µ and the angular confidence interval I ϕ are retrieved based on the position and the width of the main peak in M ϕ (using the first and second moments of the distribution).The analysis of M ϕ involves a peak-detection routine based on a multi-Gauss parameterisation.The latter is done to ensure that the retrieved parameters ϕ µ and I ϕ are not falsified due to potential additional peaks in the distribution (e.g. a cloud passing the scene, e.g.illustrated in Fig. A1 (5) The projected plume velocity vector for the ROI can then be calculated as where f and pix denote lens focal length and the pixel pitch of the detector, and d pl is the distance between the camera and the plume.Ill-constrained motion vectors in the DVF can then be identified with a certain confidence based on M ϕ and In this article, the method is demonstrated using the OpenCV (Bradski, 2000) Python implementation of the Farnebäck OF algorithm (Farnebäck, 2003, also used in Peters et al., 2015).It is pointed out, though, that it can be applied to DVFs from any motion estimation algorithm.

Applicability and uncertainties
The proposed method offers an efficient solution to identify flow vectors containing actual gas movement and separate them from unphysical results in the DVF.The method is based on a local statistical analysis of the histograms M ϕ and M |f | .A number of quality criteria were defined in order to ensure a reliable retrieval of the local displacement parameters: 1.A minimum fraction r min of all pixels in the considered ROI is required to exceed the minimum magnitude |f | min .The latter can, for instance, be set equal to one or can be estimated based on the flow vector magnitudes retrieved in a homogeneous image area (e.g.randomly oriented sky background areas in Fig. 3).
2. The same minimum fraction r min of pixels is required to match the angular expectation range specified by I ϕ (at a certain confidence level nσ ; cf.Appendices A and B3).
3. If additional peaks are detected in M ϕ , they are required to stay below a certain significance value S. The latter is measured relative to the main peak based on the integral values (cf.Appendix A3 and Fig. A1).
If any of these constraints cannot be met, the analysis is aborted.The settings used in this study are summarised in Table B1.
Please note that the method cannot account for any uncertainties intrinsic to the used OF algorithm, since these directly propagate to the derived histogram parameters.It is therefore recommended to assess the performance of the used OF algorithm independently and before applying the histogram correction (see e.g.Baker et al., 2011, Menze andGeiger, 2015).The Farnebäck algorithm used in this study showed sufficient performance both in Peters et al. (2015) and in the KITTI benchmark (cf.Menze and Geiger, 2015).The latter find that the algorithm yields correct velocity estimates in about 50 % of all cases (approximately 1σ ).Here, "correct" means that the disparity between a retrieved flow vector endpoint and its true value does not exceed a threshold of 5 %.We therefore assume that the majority (i.e.≈ 3σ ) of all successfully constrained flow vectors lie within a disparity radius of 15 %.Based on this, we assume an intrinsic, conservative uncertainty of 15 % for the effective velocities (Eq.4) retrieved from successfully constrained flow vectors.Note that this is a somewhat arbitrary choice of the intrinsic uncertainty of the Farnebäck algorithm, solely based on the findings of Menze and Geiger (2015).However, we remark again that it is beyond the scope of this paper to verify the accuracy of the Farnebäck algorithm, which we use to illustrate the performance of our new post-analysis method.For all illconstrained motion vectors which are replaced by the PDV, we assume a conservative uncertainty based on the width nσ of the histogram peaks (cf.Appendices B3 and B5.1).
Finally, we point out that the proposed histogram correction does not constitute any significant additional computational demands.For our data (i.e.1344 × 1024 pix) and on an Intel i7, 2.9 GHz machine, the required computation time for the correction is typically less than 0.1 s.In contrast, the  Farnebäck OF algorithm itself typically requires 1.5 s (same specs.)and can be considered fast in comparison with other solutions (e.g.Baker et al., 2011).

Results
The new method was applied to the Etna and Guallatiri data sets introduced in Sect.2.2.SO 2 emission rates (Eq.3) of both sources were retrieved as described in Sect.2.3 along the corresponding PCS lines (cf.Figs. 1 and 2).In order to assess the performance of the proposed correction we use the following three methods to estimate the gas velocities in the plume: 1. glob is based on a cross-correlation analysis at the position of the PCS line (i.e. the estimated velocity is applied to all pixels on and to all images of the time series).2. flow_raw uses the raw output from the Farnebäck algorithm (i.e.without correction for erroneous flow vectors).
3. flow_hybrid uses reliable optical flow vectors, identifies and replaces unphysical vectors using the DVF from the histogram analysis.
Table B1 (in Appendix B3) summarises all relevant settings for the OF-based velocity retrievals.Note that the required minimum magnitude for successfully constrained motion vectors was set per image and ROI using the lower end of I |f | at 1σ confidence.In order to assess the impact of unphysical motion vectors on the retrieved SO 2 emission rates, we define the ratio κ: where χ all corresponds to the SO 2 integrated column amount (ICA) considering all pixels on while χ pix ok corresponds to the SO 2 ICA considering only pixels showing reliable flow vectors.κ = 1, for instance, means that all motion vectors on are considered reliable.The ROIs for the OF histogram analysis were defined for each PCS line individually (based on the position and orientation of the line).

Etna results
The OF gas velocities in the plume were calculated from onband OD (τ on ) images, since the OF algorithm showed best performance for the on-band OD images (based on visual inspection before the analysis).Figure 4 shows an example DVF of the Etna plume and the corresponding histograms M ϕ and M |f | .Along the orange line, the OF algorithm performs considerably well with only 7 % of the velocity vectors found ill-constrained.If not corrected for, these unphysical motion vectors would result in an underestimation of only 1 % in the SO 2 emission rates.For the blue line, on the contrary, a total of 45 % of the pixels on were found unreliable.Moreover, many of these are located in regions showing large SO 2 -CDs.Hence, the impact is considerably large and, if not corrected for, would induce an underestimation of 33 % in the SO 2 emission rates.
Prior to the emission-rate analysis, the proposed histogram method was applied to all τ on images in order to retrieve time series of the four correction parameters p = (ϕ µ , ϕ σ , |f | µ , |f | σ ).Missing data points (i.e.where the required constraint parameters were not met; cf.Sect.2.4.1) were interpolated.The results were averaged in time using a combined median filter of width 3 (to remove outliers) and a Gaussian filter (σ = 5, to remove high-frequency variations in the retrieved DVFs).The results of this pre-analysis Figure 6 shows the results of the SO 2 emission-rate analysis for both PCS lines and the three different velocity retrieval methods.Also included are the corresponding effective velocities (average along , second panel) and the retrieved κ values (Eq.7).The latter indicates the percentage impact of unphysical OF motion vectors on the SO 2 emission rates.The plotted uncertainties in the SO 2 emission rates and the effective velocities (shaded areas) were calculated as described in Appendix B5.SO 2 emission rates between 4.9 and 9.7 kg s −1 (average of 7.1 kg s −1 ) and 4.8-10.7 kg s −1 (average of 7.8 kg s −1 ) were retrieved along the orange and blue lines, respectively, using the proposed flow_hybrid method.The slightly higher values in the aged plume are likely due to the fact that this line captures more of the emissions from the other Etna craters (cf.Supplement video no. 1).The corrected OF emission rates show good agreement with the results using the crosscorrelation velocities (glob method).The latter, however, tend to be slightly increased by about +14 % (cf.Fig. 7).The flow_raw method (i.e.uncorrected OF velocities), on the contrary, often yields significantly decreased SO 2 emission rates, especially in situations where unphysical OF motion vectors coincide with either of the retrieval lines (i.e.low κ value; cf.Fig. 4).The latter show rather strong fluctuations between consecutive frames (i.e.local scatter in the κ values) with an average impact of κ = (0.68 ± 0.15).These fluctuations are due to the somewhat random nature of the initial problem.Namely, that the occurrence (and position) of regions containing unphysical motion vectors can change significantly between consecutive frames (cf.Fig. 4).These unphysical fluctuations (in the estimated gas velocities) directly propagate to the SO 2 emission rates (retrieved using the flow_raw method) and are thus not to be misinterpreted with actual (high-frequency) variations in the SO 2 emission rates.
Relative deviations of the three methods are shown in Fig. 7 (normalised to the results from the proposed flow_hybrid method).The cross-correlation-based retrievals (glob) tend to yield slightly larger SO 2 emission rates (by +14 % on average), while the uncorrected OF (flow_raw) often shows underestimated results (by −20 % on average).However, we point out again that these underestimations generally show a rather strong variability.This includes cases showing considerably large underestimations (up to 62 %) and other cases in which the OF algorithm appears to perform sufficiently well (i.e. = 1 in Fig. 7).

Guallatiri results
The OF gas velocities for the Guallatiri data were retrieved using the on-band OD images.An example DVF is shown in The time series of the interpolated and smoothed displacement parameters for both PCS lines is shown in Fig. 9. Compared to Etna, the two plumes show considerably more variability both in orientation and in the velocity magnitudes (cf.Figs. 5 and 9).The resulting average values are ϕ µ = (12.6 ± 16.8) • and |f | µ = (1.17± 0.33) pix s −1 (crater) and ϕ µ = (15.9± 13.1) • and |f | µ = (1.30± 0.13) pix s −1 (fumaroles).Due to the rather strong temporal variations, the emissions of both sources could not always be successfully separated using the two (fixed) PCS lines.This can be seen in the Supplement video no. 2, which shows the evolution of SO 2 emission rates for both PCS lines.1. Emission rates were retrieved using the three different velocity retrieval methods described above.Uncertainties (shaded areas) are only plotted for the flow_hybrid method and the cross-correlation method (glob).Also included are time series of effective velocities (Eq.4, middle panel) and κ values (Eq.7, bottom) retrieved from the proposed histogram analysis.
Figure 7. Relative deviations of retrieved SO 2 emission rates shown in Fig. 6 for the "young_plume" (a) and the "aged_plume" (b) PCS lines using the same colour codes as in Fig. 6.The ratios are plotted relative to the results of the proposed flow_hybrid method.Results based on the cross-correlation analysis tend to be slightly larger (by about +14 %), while the uncorrected OF velocities often yield underestimated SO 2 emission rates (up to 62 %).
The results of the emission-rate analysis are shown in Fig. 10, again, including effective velocities and κ values for both PCS lines (cf.Fig. 6).As in the Etna example, the SO 2 emission rates were calculated using the three different velocity retrieval methods introduced above (i.e.glob, flow_raw, flow_hybrid).In general, similar trends can be observed.The uncorrected OF often causes significant underestimations in the SO 2 emission rates.It furthermore accom-panies rather strong (and unphysical) high-frequency fluctuations which are propagated to the SO 2 emission rates (see Sect. 3.1 for a discussion).The cross-correlation velocity analysis could only successfully be applied to the emissions from the fumarolic field (cf.Fig. 2 and Sect.B4), since the central crater plume showed too strong fluctuations both in space and time.The corresponding emission rates of the fumarole emissions (Fig. 10, right, purple) show good agreement with the flow_hybrid method.
The SO 2 emission rates, which were calculated based on the proposed flow_hybrid method, show only little variation in the central crater emissions with values ranging between 0.1 and 1.5 kg s −1 (Fig. 10, left).The corresponding fumarole emissions, however, show rather strong variations with peak emission rates of 2.5 kg s −1 (at 14:55 UTC), even exceeding the observed amounts from the central crater.The sum of both sources yields total SO 2 emission rates of tot = 1.3 ± 0.5 kg s −1 with peak emissions of up to 2.9 kg s −1 .
Relative deviations of the retrieved SO 2 emission rates between the three velocity methods are shown in Fig. 11.As in the case of Etna, the cross-correlation-based results (glob, fumaroles) tend to be slightly increased (here +23 %), while the uncorrected OF (flow_raw) results in an average underestimation of −20 %.

Corrected gas velocities
The proposed histogram correction was successfully applied to the two example data sets from Mt Etna and Guallatiri.In particular, the rather turbulent Guallatiri case clearly demon-   strated the necessity of localised gas velocity retrievals (both in the spatial and temporal domain).We showed that the Farnebäck OF algorithm is (generally) able to resolve the 2-D velocity fields in great detail.However, we also showed that unphysical OF motion vectors often induce significant underestimations in the retrieved SO 2 emission rates.For both data sets, the proposed histogram method was able to account for this issue and resulted in more robust SO 2 emission-rate retrievals (cf.Figs. 7 and 11).The corrected results show good coincidence with SO 2 emission rates based on the assumption of a global constant velocity (retrieved using a crosscorrelation algorithm).However, the limitations of the crosscorrelation method were clearly demonstrated in the case of the turbulent Guallatiri plume.

Retrieved SO 2 emission rates
The retrieved Etna emission rates of ∼ 8 kg s −1 (∼ 700 t d −1 ) are at the lower end of typically observed values (> 1000 t d −1 , e.g.Salerno et al., 2009), ranging from a few up to several hundred kg s −1 , depending on the activity (e.g.Edner et al., 1994, D'Aleo et al., 2016).The comparatively low values may be partly due to the fact that mainly the emissions of the NEC were captured, which nonetheless appeared to be the strongest source during the observation period.The presented time series of about 15 min duration is too short to infer any reliable conclusions related to the state of activity.Nonetheless, it may be noted that the measurements Figure 11.Relative deviations of Guallatiri emission rates shown in Fig. 10.The deviations are plotted as ratios normalised to the results from the proposed method (flow_hybrid), both for the crater (a, no cross-correlation results available) and for the fumarolic emissions (b).The average ratios are 1.23 ± 0.32 % (fumaroles, glob) and 0.85 ± 0.12 (crater, flow_raw) and 0.75 ± 0.22 % (fumaroles, flow_raw).Again, the latter show a rather strong variability between the images.
were recorded about 2 months prior to a major eruptive event, and that indications of decreased pre-eruptive SO 2 emissions have been observed before at Mt Etna (e.g.Caltabiano et al., 1994).A longer record of Etna's SO 2 emissions (i.e. during the months of September-December 2015) would hence be desirable in order to evaluate whether these low emission rates were characteristic of the time period prior to the eruption.
In the case of Guallatiri, these are the first SO 2 emission rates reported in the literature.This makes the retrieved emission rates of ∼ 1.5 kg s −1 (peak of ∼ 3.0 kg s −1 ) an important finding of this study.However, also in this case, the presented time series is rather short and hence not suited to inferring typical emission characteristics of the volcano.Future investigations are highly desirable in order to infer more detailed information related to the emission characteristics of Guallatiri (e.g.long-term averages).Conducting these is obviously more challenging than in the Etna case due to the remote location of the volcano.Space-based observations of this considerably weak source may be an option in the future but appear to be difficult with currently available instrumentation (e.g.Carn et al., 2016).

Conclusions
In this article, a new method for image-based gas velocity measurements in plumes was presented.The success of the method lies in the extraction of quantitative information about gas dynamics inside an emission plume by using the physical information present in a remotely recorded videoimage sequence.The method is based on local gas velocity retrievals using optical flow (OF) algorithms.OF algorithms are a powerful tool for measuring the plume velocities in great detail.However, an intrinsic weakness of such algorithms is that they rely on contrast in the images.Hence, they often yield unphysical motion estimates in low-contrast image regions (e.g. in the centre of an extended, homoge-neous plume).We showed that this weakness is unacceptable for applications relying on accurate velocity measures at specific image coordinates (such as the discussed application of camera-based SO 2 emission-rate retrievals).
The proposed method aims to address this issue based on a local post-analysis of an OF displacement vector field (DVF).The central idea is to separate reliable from unreliable motion vectors in the DVF.This is done based on distinct peaks in histograms of the DVF, allowing the local average velocity vector to be derived.The latter can then be used as a replacement for unphysical results in the DVF.The relevance of the correction was discussed using the example of SO 2 emissionrate retrievals from UV camera data.Specifically, the SO 2 emissions of Mt Etna (Italy) and Guallatiri (Chile) were investigated using two short example data sets (of about 15 min each).The gas velocities were analysed using the Farnebäck OF algorithm.Based on these data we find that unphysical motion vectors occur rather frequently and hence often induce significantly underestimated SO 2 emission rates.We further show that the correction provides an efficient solution to this problem, resulting in more reliable velocity estimates and hence in more robust SO 2 emission-rate retrievals.The proposed method therefore provides an important and useful extension for OF-based gas velocity retrievals.
1. Extraction of all displacement vectors f = [ i, j ] T within a specified ROI: F = {f | f ∀(i, j ) ∈ ROI}, where i, j denote pixel coordinates in the image.Note that the considered pixels in the ROI may be further restricted, for instance by applying an intensity threshold (e.g.we use a τ on threshold to identify plume pixels).
2. Determination of magnitude |f | = i 2 + j 2 , and orientation angle ϕ(f ) = atan2( i, j ) of all vectors in F . 5. Perform multi-peak analysis of M ϕ using a multi-Gauss fit (for details see next Sect.A2).
6. Use multi-Gauss fit result to check whether an unambiguous peak can be identified in M ϕ .If this is the case, estimate the expectation interval I ϕ = [ϕ µ − nϕ σ , ϕ µ + nϕ σ ] from the first and second moments of M ϕ (with n specifying a certain confidence level).

A2 multi-Gauss fitting routine
The multi-Gauss fitting routine is used to detect and parameterise distinct peaks in a given orientation histogram M ϕ calculated from a DVF.The parameterisation is performed by fitting a number K of superimposed Gaussians of the form to M ϕ , with the normal distribution and the corresponding parameter vector p = (p 1 , . .., p K ) = ((A 1 , µ 1 , σ 1 ), . .., (A K , µ K , σ K )).
In order to achieve a physically more reliable result in the optimisation, we recommend restricting the individual p k to certain expectation ranges: minimum required amplitude: A k > A k,min (e.g. to avoid fitting all noise peaks), lower threshold for standard deviation: (FWHM must equal or exceed histogram bin resolution), peak position in angular range: µ k ∈ {−180, 180 }, defined upper limit for allowed number of superimposed Gaussians: K ≤ K max .
A routine to perform this fit was written in Python and is implemented in the software package Pyplis (class Multi-GaussFit).The algorithm aims to find the minimum number K of Gaussians required to meet the constraint C : R pp (i) ≤ A k,min , where R pp (i) is the peak-to-peak value of the current fit residual R(i) = f K (ϕ, p) − M ϕ at iteration step i.If no additional peaks are found in R(i), the latest optimised parameter vector p is assumed sufficient.Else, p is extended by all additionally detected peaks (in R(i)) and the leastsquares fit is re-applied until the optimisation constraint C is fulfilled or until a break constraint is met (e.g.maximum iteration reached or maximum number of allowed Gaussians).An exemplary fit result is shown in Fig. A1.
A3 Retrieval of main peak parameters from multi-Gauss fit result The multi-Gauss parameterisation of M ϕ allows us to identify the most prominent peak in the distribution (which may be a superposition of several Normals) and separate it from potential additional peaks.The latter can have significant impacts on the retrieved statistical parameters ϕ µ , ϕ σ (cf.Fig. A1).
Numerically, the retrieval of the main peak parameters from a given fit result vector p is performed as follows: 1.For a given Gaussian p k in p, find all fitted Gaussians within a specified confidence interval (nσ ) k around p k and calculate the integral value I k of the local overlap (I k corresponds to the number of vectors belonging to the main peak).
2. Do step 1. for all detected Gaussians in p, resulting in a vector I K of length K containing integral values of the local overlaps.
3. Find the main peak position based on the index k * showing the largest integral value (in I K ).
5. Retrieve the mean and standard deviation of p based on the first and second moments of the resulting main peak distribution f K (ϕ; p ) (cf.Eq.A1), i.e.
resulting in an estimate of the predominant displacement direction in the ROI: Appendix B: Emission-rate analysis supplementary information B1 DOAS SO 2 retrieval The DOAS spectra from both data sets (Etna, Guallatiri) were analysed using the software DOASIS (Kraus, 2006).All spectra were corrected for electronic offset and dark current and were analysed using a clear sky Fraunhofer reference spectrum (FRS) recorded close in time (to keep potential O 3 interferences at a minimum).In addition, a Ring spectrum, determined from the FRS, was fitted as well as the absorption cross sections (XS) of SO 2 (Hermans et al., 2009) and O 3 (Burrows et al., 1999).The latter were convolved with the instrumental line spread function (using the measured 334.15 nm mercury line).FRS and Ring were linked to each other and were allowed a slight shift of 0.2 nm and squeeze of 2 %.The same shift and squeeze was allowed for the two XS, which were also linked.The retrieval was performed between 314.6 and 326.4 nm.A third-order DOAS polynomial was fitted to account for broadband extinction and an additional zero-order offset polynomial (fitted in intensity space) was included to account for instrumental effects (e.g.stray light).

B2 Camera calibration
Figure B1 shows the DOAS calibration curves retrieved for both data sets.The camera AA values correspond to the pixels covered by the DOAS-FOV shown in Figs. 1 and 2. Prior to the calibration, the camera images and the DOAS data were merged in time.Note that the calibration data displayed in Fig. B1 is not corrected for the signal-dilution effect.In order to calibrate the dilution-corrected AA images, the fitted calibration polynomial was extrapolated into the AA regime of the dilution-corrected images.This is based on the assumption that the calibration curve also remains linear at larger optical densities, which is justified by the considerably good plume conditions (low to no condensation) and the low to moderate range of observed SO 2 -CDs (cf.Sect.2.3.3;see also Kern et al., 2013).Furthermore, this calibration method assumes that the retrieved DOAS SO 2 -CDs exhibit approximately the same amount of signal dilution as the camera imagery.This is justified by the fact that the DOAS analysis was applied in a wavelength range coinciding with the on/offband regime of the camera filters (cf.previous Sect.B1).

B3 Settings for optical flow retrieval
All relevant settings for the optical-flow-based gas velocity retrievals are summarised in Table B1

B4.2 Guallatiri
Figure B3 shows the result of the velocity cross-correlation analysis using the blue (fumarole) PCS retrieval line (cf.Sect.3.2) resulting in a gas velocity of v glob = 3.49 m s −1 .

B5 SO 2 emission-rate uncertainties
Uncertainties in the presented emission rates (shaded areas in Figs. 6 and 10 top) were calculated based on Eq. (3) using Gaussian error propagation.Uncertainties in the plume distance (from uncertainty in plume azimuth and camera viewing direction), in the retrieved SO 2 -CDs (from slope error the calibration polynomial) and in the effective gas velocities (Eq.4) were considered.The latter were assumed constant for cross-correlation-based velocities using v glob = 1.5 m s −1 .For the optical-flow-based retrievals, the uncertainties were estimated per image and PCS line as described in Sect.2.4.1.Note that uncertainties resulting from potential radiative transfer effects were not included.These are discussed in Sect.2.3.3.Corresponding time series of integrated on-band ODs (dashed blue and magenta), and further, the shifted PCS signal corresponding to the retrieved correlation lag of 18.0 s (solid blue).Note that here, the velocity analysis was applied using a time series of on-band OD images rather than the τ AA which was used in Fig. B2.

B5.1 Sensitivity to the chosen histogram analysis settings
The sensitivity of the retrieved emission rates (cf.Sect.3) to the input parameters |f | min and nσ of the proposed histogram correction method (cf.Appendix B3) was investigated.These two parameters have the largest impact on the results since they determine which flow vectors are considered ill-constrained and which ones are not.The sensitivity analysis was performed using the proposed flow_hybrid method applied to 30 images from the Etna data set that were not corrected for the signal-dilution effect, since the latter is irrelevant for this study (all other analysis settings are the same as described in Sect.2.3). Figure B4 shows the results of this investigation.The choice of n has a rather small impact on the emission rates, whereas the choice of |f | min impacts within a range of approximately ± 17 %.However, considering the more realistic interval of 1-2 for |f | min , the impact is less than 8 %.The investigated value ranges are 0-4 pixels for |f | min and 1-3 for the n and the deviations are plotted as percentage deviations SO 2 from the average SO 2 emission rate retrieved from this study.The latter amounts to 2.4 kg s −1 (not dilution corrected; see text).The analysis was performed using the average SO 2 emission rates retrieved from 30 images of the Etna data set.

Figure 1 .
Figure 1.(a) Etna overview map showing position and viewing direction of the camera (camera cfov, fov) which was located on a rooftop in the town of Milo.Also indicated is the summit area (source) and the plume azimuth (plume direction).(b) Example SO 2 -CD image of the Etna plume including two PCS lines (orange, blue) used for emission-rate retrievals and two corresponding offset lines (green, red) that are used for cross-correlation-based plume velocity retrievals (cf.Appendix B4).Position and extent of the DOAS-FOV for the camera calibration are indicated by a green spot.Note that the displayed plume image is size reduced by a factor of 2 (Gauss pyramid level 1).

Figure 2 .
Figure 2. (a) Guallatiri overview map showing position and viewing direction of camera (camera cfov, fov), summit area (source) as well as the plume azimuth (plume direction).(b) Example SO 2 -CD image of the Guallatiri emissions including two PCS lines used to retrieve SO 2 emission rates from the central crater (orange) and from a fumarolic field (blue) located behind the flank in the viewing direction.An additional line (magenta) is used to estimate gas velocities using a cross-correlation algorithm (relative to blue PCS line; cf.Appendix B4).The position and extent of the DOAS-FOV are indicated by a green spot.Note that the displayed plume image is size reduced by a factor of 2 (Gauss pyramid level 1).
Figure 3 shows an example DVF (left) retrieved from the Etna plume including an example rectangular ROI (top).Two further images show the corresponding OF displacement orientation angles ϕ (middle) and flow vector magnitudes |f | (bottom).Histograms M (i.e.M ϕ , M |f | ) of the motion field are plotted in the right panels and were calculated considering all image pixels belonging to the displayed ROI.From the images and histograms, certain characteristics become clear: 1. Image regions containing unphysical motion estimates are characterised by (local) random orientation and short flow vectors (cf.sky background pixels).2. These unphysical motion vectors manifest as a constant offset in M ϕ and as a peak at the lower end of M |f | .3. Image regions showing reliable motion estimates, on the other hand, are characterised by (locally) homogeneous orientation ϕ and magnitudes |f | exceeding a certain minimum length |f | min .4. These successfully constrained motion vectors manifest as distinct peaks in M ϕ and M |f | .5.The width of these peaks can be considered a measure of the local fluctuations or the variance of the velocities (e.g. a very narrow and distinct peak in M ϕ would indicate a highly directional movement).
).Based on the analysis of M ϕ , a second histogram M |f | is determined, containing the displacement magnitudes |f | of all vectors matching the angular confidence interval I ϕ and exceeding the required minimum magnitude |f | min .Also here, an expectation value |f | µ and confidence interval I |f | are estimated based on the first and second moments of the histogram.The analysis yields four parameters p ROI = (ϕ µ , ϕ σ , |f | µ , |f | σ ) which are used to calculate the PDV within the corresponding ROI:

Figure 3 .
Figure 3. (a) Example output of the Farnebäck optical flow algorithm including a rectangular ROI.Two further images show the corresponding orientation angles ϕ (b) and magnitudes |f | (c) of the DVF.Corresponding histograms M ϕ and M |f | are plotted on the right and include all pixels in the displayed ROI.The histograms are plotted both including (dashed lines) and excluding (red and blue shaded areas) short flow vectors (i.e.|f | > |f | min = 1.5pix.The orientation angles are plotted in an interval −180 • ≤ ϕ ≤ +180 • where −90 • , +90 • correspond to −i, +i directions and 0 • to the vertical upwards direction (−j ).The DVF was calculated using two consecutive AA images ( t = 4.0s) of the Etna plume, recorded on 16 September 2015 at 07:14.

Figure 4 .
Figure 4. (a) Example flow vector field (blue lines with red dots) of the Farnebäck optical flow algorithm for the Etna plume at 07:13 UTC including the two PCS lines (blue and orange) and the corresponding ROIs used for the histogram analysis (semi-transparent rectangles).Middle, right: histograms of orientation angles M ϕ and vector magnitudes M |f | for both lines (bar plot), determined using condition no.7 in Appendix A1.The M ϕ histogram (b) also includes fit results of the multi-Gauss peak detection (thick solid lines).The retrieved histogram parameters (ϕ µ , |f | µ ) and expectation intervals M ϕ , M |f | are indicated with solid and dashed vertical lines, respectively.From the corrected DVF, average effective velocities of v eff = (3.9± 0.5) m s −1 (orange line) and v eff = (4.4± 0.8) m s −1 (blue line) were retrieved.Note that in the left image (1) vectors shorter than 1.5 pixels are excluded, (2) the displayed vector lengths were extended by factors of 3, and (3) only every 15th pixel of the DVF is displayed.

Figure 5 .
Figure 5.Time series of retrieved PDV parameters ϕ µ and |f | µ (dashed lines) for the two Etna PCS lines (same colours; cf.Fig. 1) and corresponding values after applying interpolation and smoothing (solid lines).The expectation intervals I ϕ and I |f | are plotted as shaded areas.

Fig. 8 .
Here, the two sources are clearly separable, showing a convective central crater plume (approx.location at cols.i ≈ 50-80) and the emissions from the fumarolic field located behind the volcanic flank (i ≈ 100-300).Also included are the results of the proposed OF histogram analysis, which was performed relative to the two displayed PCS lines used for the SO 2 emission-rate analysis (cf.Fig.2).In this example, the OF algorithm performed considerably well.The uncorrected OF would therefore result in small underestimations of 6 % (crater) and 3 % (fumaroles) in the SO 2 emission rates.The different plume characteristics of both sources can be clearly identified based on the displayed histogram distributions.The central crater plume (orange) rises almost vertically (ϕ = (−12.4± 17.5) • ) and reaches velocity magnitudes of up |v| max = 4.5 m s −1 .The fumarolic emissions (blue) are less convective (ϕ = (+55.6± 17.4) • ) and show slightly smaller velocities with maximum magnitudes of |v| max = 3.4 m s −1 .

Figure 6 .
Figure 6.Time series of Etna emission rates (top panel), showing emissions of the young (a) and the aged (b) plume of the NE crater (orange and blue) using the two PCS lines shown in Fig.1.Emission rates were retrieved using the three different velocity retrieval methods described above.Uncertainties (shaded areas) are only plotted for the flow_hybrid method and the cross-correlation method (glob).Also included are time series of effective velocities (Eq.4, middle panel) and κ values (Eq.7, bottom) retrieved from the proposed histogram analysis.

Figure 8 .
Figure 8.(a) Example output of the Farnebäck optical flow algorithm for the Guallatiri emissions at 14:48 UTC including the two example PCS lines (blue/orange line) and the corresponding ROIs.(b, c) Histograms of magnitudes M ϕ and orientation angles M |f | used to retrieve the expectation intervals I ϕ and I |f | and the corresponding PDV in each ROI (cf.Eq. 5).From the latter, effective velocities of v eff = (3.1 ± 0.5) m s −1 (crater, orange) and v eff = (1.8 ± 0.6) m s −1 (fumaroles, blue) were retrieved.Note that in the left image (1) vectors shorter than 1.5 pixels are excluded, (2) the displayed vector lengths were extended by a factor of 2, and (3) and only every 15th pixel of the DVF is displayed.

Figure 9 .
Figure 9.Time series of retrieved PDV parameters ϕ µ and |f | µ (dashed lines) for the two Guallatiri PCS lines (same colours; cf.Fig. 8) and the corresponding values after applying interpolation and smoothing (solid lines).The expectation intervals I ϕ and I |f | are plotted as shaded areas.

Figure 10 .
Figure10.Guallatiri SO 2 emission rates from the summit crater (a, orange) and the fumarolic field (b, blue) using the two PCS lines shown in Fig.2.Uncertainties (shaded areas) are only plotted for the flow_hybrid and the glob (cross-correlation-based) velocity-retrieval methods.Also included are the corresponding effective velocities (from the flow_hybrid method) and the OF quality factors κ (Eq.7).The central crater emissions show only little variability ( ≈ 0.6 kg s −1 ), while the fumarolic emissions are characterised by a comparatively strong emission event at 14:55 UTC showing peak emissions of 2.5 kg s −1 .
3. Extraction of all vectors in F exceeding a certain magnitude threshold |f | min : F = {f : f ∈ F ∧ |f | > |f | min }. 4. Calculation of orientation angle histogram M ϕ (F , ϕ) of vectors in F (with the histogram bin-width ϕ).

7.
Extraction of all flow vectors matching the angular expectation intervalI ϕ : F = {f : f ∈ F ∧ ϕ(f ) ∈ I ϕ }.8.Calculation of a displacement length histogram M |f | (F , |f |) from vectors in F (with |f | being the bin-width in units of pixel displacements).9. Determine average displacement length |f | µ and standard deviation |f | σ using first and second moments of M |f | .

Figure
Figure A1.(a) Example fit result (blue line) of a multi-peak analysis applied to synthetic data (blue crosses).(b) Corresponding fit residual.The data set consists of three Gaussian Normals including Gaussian noise (with σ noise = 9 counts).Two overlapping peaks are located at µ = −110 (A = 150, σ = 25) and µ = −50 (A = 300, σ = 20) (forming the predominant peak) and one additional peak at µ = 90 (A = 150, σ = 10) with a significance of 16 %.Expectation parameters (ϕ µ , ϕ σ ) are indicated with solid and dashed vertical lines, respectively.The latter were retrieved as described in Appendix A3, both including and excluding the additionally detected peak at µ = 90, in red and green.

Figure B1 .
Figure B1.DOAS calibration curves of Etna data (a) and Guallatiri data (b) retrieved from AA images (not dilution corrected) within the image region covered by the DOAS-FOV (cf.Figs. 1 and 2).The corresponding SO 2 -CDs retrieved with the DOAS instrument are plotted on the y axis.Y axis offsets were corrected for during the calibration of the AA images.

Figure B2 .
Figure B2.Result of velocity cross-correlation analysis for the Etna data.(a) Example plume AA image including two PCS lines (orange and blue) and corresponding offset lines (green and red) used for the analysis.(b, c) Corresponding time series of the integrated AA values along the two PCS lines (original: dashed, shifted using correlation lag: solid) and along the offset lines (solid) using the same colour scheme.

Figure B3 .
Figure B3.Result of velocity cross-correlation analysis for the Guallatiri data.(a) Example plume on-band OD image (τ on ) including the PCS (blue) and offset line (magenta) used for the analysis.(b)Corresponding time series of integrated on-band ODs (dashed blue and magenta), and further, the shifted PCS signal corresponding to the retrieved correlation lag of 18.0 s (solid blue).Note that here, the velocity analysis was applied using a time series of on-band OD images rather than the τ AA which was used in Fig.B2.

Figure B4 .
Figure B4.Sensitivity of SO 2 emission rates as a function of the chosen input settings |f | min (y axis) and the confidence level n (x axis).The investigated value ranges are 0-4 pixels for |f | min and 1-3 for the n and the deviations are plotted as percentage deviationsSO 2 from the average SO 2 emission rate retrieved from this study.The latter amounts to 2.4 kg s −1 (not dilution corrected; see text).The analysis was performed using the average SO 2 emission rates retrieved from 30 images of the Etna data set.

Table 1 .
Instrumental set-up during both campaigns.