A novel simplified ground-based thermal infrared (TIR) system for volcanic plume geometry, SO₂ columnar abundance, and flux retrievals

In the last few decades, volcanic monitoring using remote sensing systems has become an essential tool to investigate the effects of volcanic activity on environment, climate, human health, and aviation, as well as to give insights into volcanic processes. Compared to satellite measurements, ground-based instruments offer continuous spatial and temporal coverage capable of providing high-resolution and high-sensitivity data.

This work presents a new simplified prototype of a thermal infrared (TIR) system (named “VIRSO2”). The instrument comprises three cameras, one working in the visible spectrum and two in the TIR (8–14 µm). In front of one of the two TIR cameras, an 8.7 µm filter is placed. The system is designed for the detection of volcanic emission, geometry estimation, determining the columnar content of SO₂ and ash, and SO₂ flux retrievals. The retrieval procedures developed are detailed starting from the geometric characterisation with wind direction correction, the calibration by considering the effects of filter multireflections and temperature, and the SO₂ mass by exploiting MODTRAN radiative transfer model (RTM) simulations. The SO₂ flux is then computed by applying the traverse method, with the plume speed obtained from the wind speed at the crater altitude. As test cases, the measurements collected at Etna Volcano (Italy) on 1 April 2021 during a lava fountain episode and 30 August 2024 during a quiescent phase have been considered. The results show that the system can provide reliable information on plume detection, altitude, and SO₂ flux.

Its simplicity, its low cost, and the possibility of carrying out measurements at a safe distance from the vent during both the day and the night make this system ideal for real-time monitoring of volcanic emissions, thus helping to provide information on the state of activity of the volcano and therefore to mitigate the effect that these natural phenomena have on humans and the environment.

1 Introduction

During their activities, volcanoes emit large quantities of particles (ash, water droplets) and gases (mostly H₂O, CO₂, SO₂, and HCl) into the atmosphere, with severe threats to the environment (Thordarson and Self, 2003; Craig et al., 2016), climate (Robock, 2000; Haywood and Boucher, 2000; Grainger and Highwood, 2003; Solomon et al., 2011; von Savigny et al., 2020), human health (Horwell and Baxter, 2006; Horwell et al., 2013), and aviation (Casadevall, 1994; Kristiansen et al., 2015). In the last few decades, ground and satellite remote sensing systems have become essential tools for monitoring volcanic activity in real time to mitigate the effect of these natural phenomena and give insight into volcanic processes. If satellite measurements guarantee global spatial coverage with generally low spatial resolution and variable accuracy, ground systems offer continuous spatial and temporal coverage capable of providing high resolution and sensitivity but are limited to specific locations.

From satellite, the volcanic clouds are generally detected in the thermal infrared (TIR) spectral range by exploiting the opposite spectral behaviour between ash and water particles at around 11 and 12 µm (Prata, 1989a; Prata, 1989b; Wen and Rose, 1994) and the wide SO₂ absorption features at 7.3 and 8.7 µm (Realmuto et al., 1994; Watson et al., 2004). The absorption feature at 7.3 µm is strongly affected by the atmospheric water vapour, so it is generally used when a volcanic cloud is located at high altitudes, while 8.7 µm lies in an atmospheric window and then can be exploited for any volcanic cloud altitude (Watson et al., 2004; Corradini et al., 2009). Moreover, remote sensing in the TIR guarantees the possibility of carrying out monitoring during both the day and the night.

As emphasised by Prata et al. (2024), TIR ground-based remote sensing applied to studies on volcanic plumes is a relatively new application and experienced a strong increase in use following the wide availability of uncooled thermal cameras. These systems have been used to characterise volcanic activity at Stromboli (Italy) (Ripepe et al., 2004; Patrick et al., 2007; Patrick, 2007), to investigate unsteadiness in low-energy eruptions of Sabancaya Volcano (Peru) (Rowell et al., 2023), to study the dynamics of volcanic plumes at Santiaguito Volcano, Guatemala (Sahetapy-Engel and Harris, 2009), to track volcanic explosions at Reventador Volcano, Ecuador (Vásconez et al., 2022), and as part of a continuous, real-time volcanic monitoring system (Mereu et al., 2023).

In this work, we will describe a system that represents an evolution of the TIR instrument described in Prata et al. (2024). In this case, two TIR and one visible (VIS) cameras are present, and in front of one of the two twin TIR cameras, an 8.7 µm filter is placed. The data collected and properly calibrated are used to detect volcanic plumes, estimate their altitude and thickness, and estimate SO₂ slant column densities (SCDs). The SO₂ flux is also obtained by knowing the viewing geometry and by exploiting the plume speed computed from the European Centre for Medium-Range Weather Forecasts (ECMWF) Reanalysis v5 (ERA5) hourly dataset (Hersbach et al., 2023). As test cases, the measurements collected on Etna (Sicily, Italy) in the period March–April 2021 and August 2024 have been considered. With its current height of 3403 m a.s.l. (September 2024), Etna is the largest European volcano. Its intense and continuous activity (e.g. 66 lava fountains from late 2020 to early 2022) (Guerrieri et al., 2023) makes it the most active volcano in Europe and one of the most active in the world.

The paper is organised as follows. In Sect. 2 the characteristics of the TIR ground-based system are described. Sections 3, 4, and 5 present the algorithms for the camera view geometrical computation, the calibration, and the SO₂ retrieval procedures. The SO₂ flux computation with a validation and sensitivity study is described in Sects. 6 and 7, while the conclusions are presented in Sect. 8.

2 System description

A novel TIR ground-based camera system (called “VIRSO2”) has been developed in the sphere of the ESA-VISTA project (https://eo4society.esa.int/projects/vista/, last access: 7 January 2025) by AIRES Pty Ltd company (https://www.aires.space/, last access: 7 January 2025). The system consists of three cameras (see Fig. 1a): one VIS, one TIR broadband (7.8–14 µm, BB – broadband), and another similar TIR broadband with a narrowband filter in front (centred at 8.7 µm, NB – narrowband) to exploit the SO₂ strong absorption signature. Figure 1b shows a table with its main characteristics: the two TIR broadband cameras consist of a 320 × 240 uncooled microbolometer detector array (manufactured by SEEK Thermal, MOSAIC Core C3 Series, https://www.thermal.com/mosaic-core-320x240-4mm.html, last access: 7 January 2025). Figure 1c shows the spectral response functions (SRFs) of BB (blue line) and NB (red line) respectively used in the MODerate resolution atmospheric TRANsmission (MODTRAN) 5.3 (Berk et al., 2005) simulations needed for quantitative retrievals. As the figure shows, a simple top-hat function between 1282 cm⁻¹ (7.8 µm) and 714 cm⁻¹ (14 µm) was considered for BB since there is no specific information about it from SEEK Thermal. For NB, the manufacturer-supplied spectral transmittance (normalised to 1) of the 8.7 µm filter was used. The MODTRAN spectral radiances (in mW m⁻² sr⁻¹ cm⁻¹ units and at 1 cm⁻¹ steps) were weighted by the two SRFs and then converted into brightness temperatures by inverting the Planck function considering a central wavelength of 998 cm⁻¹ (10.02 µm) and 1151 cm⁻¹ (8.69 µm) respectively.

Figure 1 (a) View of the system with detail of the VIS and TIR cameras and 8.7 µm filter. (b) Specifications for the SEEK Thermal IR camera (MOSAIC Core C3 Series), taken from Prata et al. (2024). (c) Spectral response functions used for broadband (BB) and narrowband (NB) channels with the SO₂ transmittance spectrum.

The system is controlled by a single-board computer with 500 GB memory for data storage, complete with Wi-Fi, Bluetooth, dual ethernet ports, USB, HDMI, and a microSD card slot to ensure adequate connectivity under most conditions. The system is operated using an external power bank. It is extremely portable (total weight of about 3 kg) and relatively low cost (about EUR 9000). Its ability to be used during both the day and the night and its good sensitivity make it a very important tool for the continuous and real-time monitoring of volcanic emissions from the ground. More information about the system can be found in Prata et al. (2024). Figure 2 shows a schematic flowchart of the procedure for SO₂ retrieval (a) and an example of a BB temperature image collected at Etna Volcano on 31 March 2021, during an eruptive event (b). The maximum value of the colour map is set to 300 K to make more evident the volcanic plume, but actually the maximum value detected by the camera is 443 K. The minimum value, relating to clear sky, is 233 K.

Figure 2 (a) Flowchart of the developed procedure for SO₂ retrieval. (b) VIRSO2 TIR broadband (BB) image collected at Etna (top of the “Sapienza” cableway, 2500 m a.s.l.) on 31 March 2021 at 21:00:04 UTC.

3 Camera geometry computation

Knowledge of the camera pixel dimensions is fundamental for the plume geometry (altitude and thickness) retrievals and SO₂ total mass or flux computations. Pixel size depends on both the focal length and the field of view (FOV) of the detector, as well as on the distance between the system and the target and the inclination of the system itself. The procedure is based on geometrical considerations with the main assumption that the focal plane is perpendicular to the ground. This approximation is generally adopted for analysis of the ascending volcanic plume in image data (Valade et al., 2014; Bombrun et al., 2018; Simionato et al., 2022). The wind orientation is an important factor to take into account in case of volcanic plumes, so it was implemented in the procedure (see Sect. 3.1).

As summarised in Fig. 1b, the TIR cameras are composed of 320 columns (NC) and 240 rows (NR), with a horizontal and vertical field of view (HFOV and VFOV) of 56 and 42° respectively. Figure 3a shows the camera measurement configuration considered. Assuming the elevation angle of the camera (α) refers to the centre of the image, it is possible to obtain the angular position (θ, φ) of the different grid points:

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathit{\theta }\left(i\right)=\mathit{\alpha }+\left({\displaystyle \frac{\text{NR}}{\mathrm{2}}}+\mathrm{1}-i\right)\cdot {\displaystyle \frac{\text{VFOV}}{\text{NR}}},\text{(2)}& {\displaystyle }\mathit{\phi }\left(j\right)=\left(j-{\displaystyle \frac{\text{NC}}{\mathrm{2}}}-\mathrm{1}\right)\cdot {\displaystyle \frac{\text{HFOV}}{\text{NC}}},\end{array}$$

where i=1, NR + 1 (from top to bottom) and j=1, NC + 1 (from left to right). In this way, θ(i) and φ(j) range between (α + 21°, α − 21°) and (−28°, +28°) respectively.

Figure 3 (a) Schematic geometrical configuration of the VIRSO2 camera FOV. The camera is placed at C. (b) VIRSO2 camera at Stromboli (Italy). (Note that Earth curvature effects are ignored as generally the system operates on small scales of ∼ 10 km or so).

The horizontal distance from the central column (x) and the vertical distance from the ground (y) of each grid point can be obtained considering D₀ (perpendicular distance between the camera and the image plane):

$$\begin{array}{}\text{(3)}& {\displaystyle }x\left(j\right)=D_{\mathrm{0}}\cdot \tan \left(\mathit{\phi }\left(j\right)\right),\text{(4)}& {\displaystyle }y\left(i\right)=D_{\mathrm{0}}\cdot \tan \left(\mathit{\theta }\left(i\right)\right).\end{array}$$

From the differences between adjacent grid points, the horizontal (dx) and vertical (dy) size, the area (a), and the height a.s.l. (h) of each pixel can be obtained (i=1, NR; j = 1, NC):

$$\begin{array}{}\text{(5)}& {\displaystyle }dx\left(j\right)=x\left(j+\mathrm{1}\right)-x\left(j\right),\text{(6)}& {\displaystyle }dy\left(i\right)=y\left(i\right)-y\left(i+\mathrm{1}\right),\text{(7)}& {\displaystyle }a\left(i,j\right)=dx\left(j\right)\cdot dy\left(i\right),\text{(8)}& {\displaystyle }h\left(i\right)=y\left(i\right)-dy\left(i\right)/\mathrm{2}+h_{\mathrm{0}},\end{array}$$

where h₀ is the altitude of the location where measurements are performed.

Figure 4 shows the horizontal (x size, Fig. 4a) and vertical (y size, Fig. 4b) dimensions of the pixels, the pixel area (Fig. 4c), and the height above sea level (Fig. 4d), computed for the “Etna-Piano del Vescovo” measurements (1 April 2021, Fig. 5).

Figure 4 Pixel dimensions and height obtained with D₀ = 6.4 km, α = 30°, h₀ = 1380 m a.s.l. (setting for the camera placed in Piano del Vescovo), and no wind correction. (a) Horizontal pixel dimension dx (m). (b) Vertical pixel dimension dy (m). (c) Pixel area (m²). (d) Height above sea level (km).

3.1 Correction for wind direction

The maps obtained in Fig. 4 are correct only if the target is exactly perpendicular to D₀, i.e. if the wind direction is parallel to the camera's focal plane. If instead the plume comes towards or goes away from the camera, it is necessary to consider that the horizontal distance from the target changes according to the wind direction and the distance from the crater position in the camera image (Fig. 5).

Figure 5 Bird's-eye schematic view representation to take into account the wind direction (WD, red line) relative to the focal plane of VIRSO2 camera.

Knowing the column position of the crater (j_C) and the relative azimuth angle (−90° < ω < 90°) between the wind direction (red line) and the focal plane of the camera (black line), it is possible to obtain, for each column of the image, the point of intersection between the lines of sight (LSs) of the camera (green lines) and the wind direction (WD) line. In a Cartesian reference system in which the camera is placed at the origin and the focal plane is horizontal, the WD and LSs are identified by

$$\begin{array}{}\text{(9)}& {\displaystyle }{m}_{\mathrm{WD}}=\tan \left(\mathit{\omega }\right),\text{(10)}& {\displaystyle }{q}_{\mathrm{WD}}=D_{\mathrm{0}}-m_{\mathrm{WD}}\cdot x\left(j_{\mathrm{C}}\right),\text{(11)}& {\displaystyle }{m}_{\mathrm{LS}}\left(j\right)={\displaystyle \frac{D_{\mathrm{0}}}{x\left(j\right)}},\text{(12)}& {\displaystyle }{q}_{\mathrm{LS}}=\mathrm{0},\end{array}$$

where m_WD and m_LS/q_WD and q_LS are the slope/axis intercept of the lines considered.

Finally, the y position of each of these intersection points gives the new $D_{\mathrm{0}}^{*}\left(j\right)$ (cyan lines in Fig. 5) to be used in Eqs. (3)–(8) to compute pixel size dimensions and heights.

Figure 6 shows two BB images of an Etna volcanic plume collected during the paroxysm of 1 April 2021 with and without the correction for the wind direction. In the upper panels, the eruption was seen during the night from “Piano del Vescovo” (37.70° N, 15.05° E; 1380 m a.s.l.; about 7 km south-east of the crater), while the bottom panels show the volcanic plume some hours later in the morning from Nicolosi (37.61° N, 15.02° E; 820 m a.s.l.; about 15 km south of the crater). Wind data were taken from the mesoscale model of the hydrometeorological service of Agenzia Regionale per la Protezione Ambientale (ARPA) Emilia Romagna, which is named ARPASIM (Scollo et al., 2009) and considers an hourly model output from a 72 h weather forecast provided for Etna every 3 h. The left panels show the plume heights obtained from Eq. (8) without wind correction (ω = 0°), while in the right panels, the effect of wind correction is considered. In the “Piano del Vescovo” example, ω = 26° and the top plume height changes from about 8.5 to 7 km. In the “Nicolosi” case, where ω = 52°, the improvements are much more important: the top plume height is now about 9–10 km a.s.l. compared to 15–16 km without wind correction. Using the Spinning Enhanced Visible and InfraRed Imager (SEVIRI) instrument, on board the Meteosat Second Generation (MSG) geostationary satellite, a volcanic plume top altitude of 6.9 km a.s.l. was estimated at 00:55 UTC (Guerrieri et al., 2023). This is in good agreement in particular with the “Piano del Vescovo” VIRSO2 nighttime measurements. In the morning, an increase in volcanic activity produced a higher volcanic plume, confirmed by INGV bulletin no. 258 at 07:42 UTC (https://www.ct.ingv.it/Dati/informative/vulcanico/ComunicatoETNA20210401074037.pdf, last access: 7 January 2025), in which it is said that the eruptive cloud exceeded 9 km a.s.l.

Figure 6 (a, b) Etna-Piano del Vescovo (1 April 2021 at 02:29:41 UTC); (c, d) Etna-Nicolosi (1 April 2021 at 07:02:04 UTC). The yellow contour lines (in km) represent the VIRSO2 pixel heights obtained without (a, c) and with (b, d) the wind direction correction.

Table 1 shows the mean (maximum) relative errors of the pixel heights varying α (from 20 to 50°), D₀ (from 2.5 to 10 km), and ω (from 0 to 60°) compared to the calibration tool developed by Snee et al. (2023). The errors increase slightly with the increase in α and D₀, while they increase more rapidly for |ω| > 50°. For wind-focal plane angles |ω| < 45°, the relative errors of all the 320 × 240 pixels of VIRSO2 are always below 1 %. So, during field measurements, it would be desirable to position the VIRSO2 so that the camera's focal plane is parallel to the wind direction or tilted at an azimuth angle of no more than 45°.

Table 1 Mean (maximum) relative errors of the pixel heights using the method presented in this paper and the calibration tool of Snee et al. (2023). The altitude of the camera is fixed at 2000 m a.s.l., and the crater is positioned in the central column of the VIRSO2 image (i_C = 200, j_C = 160).

4 VIRSO2 calibration

The calibration process is an essential procedure to obtain reliable quantitative results. Here different effects must be taken into account: at first, the non-perfect transmissivity of the 8.7 µm filter produces a “ghost image”. Then, an adjustment is also important for the BB camera, considering that the clear-sky temperature often does not match that of the MODTRAN simulations. Finally, the filter temperature strongly affects the NB measurements, so a calibration of this band is necessary.

Among the corrections said above, some minor preprocessing tasks are also needed. Firstly, a “mountain mask” is created to distinguish between land and sky pixels. This can be done manually once and for all if the camera is kept in the same position during the measurement session. Note that two separate mountain masks, one for BB and one for NB, are necessary because the two sensors, even if very close to each other, do not frame exactly the same scene. This is a common issue known as “X/Y shift” between the two cameras. Given that in the calibration procedure we have to compare BB and NB pixels in the same scene position, the two masks were also used to co-register NB to BB. Finally, a “sky mask” is created too, complementary to the “mountain mask” but also excluding a layer of 15 pixels above the ground to avoid problems with the temperature near the contour of the mountain.

4.1 Narrowband ghost image removal

All wide-field imagery contains a certain amount of optical vignetting, which is manifested by the gradual dimming of light towards the edges of the image caused by optical components obstructing off-axis light (Li and Zhu, 2009). For the BB camera, this effect is small enough for a correction to not be required (Prata et al., 2024). In addition, a so-called “ghost image” occurs in the NB channel due to the filter placement. It occurs because the interference filter is not 100 % transmissive, which results in a reflection of any hot components behind the filter back onto the detector. Since the reflected rays pass back through the focussing optics, an image can be achieved in the centre of the detector.

The correction for this effect is obtained using a purpose-built black target placed in front of the camera. The blackbody plate is large enough to fill the entire FOV, and a uniform temperature across the plate is considered due to the high proximity of the sensor to the target. The actual temperature of the plate is not required as the correction is only of a “geometric” type.

Figure 7a and b show an example of BB and NB measurements collected at “Etna-Montagnola” (37.72° N, 15.00° E; 2600 m a.s.l.) on 30 August 2024, 05:22 UTC, showing the degassing plume of Mt. Etna. The NB image (Fig. 7b) clearly shows the presence of the “ghost” effect approximately in the centre of the image.

Figure 7 TIR measurements collected at “Etna-Montagnola” on the morning of 30 August 2024. (a) BB and (b) NB original temperatures. (c) The temperature correction image that was obtained from black-target measurements. (d) NB temperatures after the removal of the ghost image.

A similar “ghost image” is present in the black-target measurement too, so a “temperature correction image” (Fig. 7c) is simply obtained as the difference between the black-target temperature image T_b,NB(ij) and its mean value $\stackrel{\mathrm{‾}}{T_{\mathrm{b},\mathrm{NB}}}$ over the 320 × 240 array:

$$\begin{array}{}\text{(13)}& {\mathrm{tci}}_{\mathrm{NB}}\left(ij\right)=T_{\mathrm{b},\mathrm{NB}}\left(ij\right)-\stackrel{\mathrm{‾}}{T_{\mathrm{b},\mathrm{NB}}}.\end{array}$$

To remove the “ghost” effect for any image, we deduct the correction factor (Eq. 13) from the temperature originally measured in each pixel:

$$\begin{array}{}\text{(14)}& T_{\mathrm{NB}}^{\prime }\left(ij\right)=T_{\mathrm{NB}}\left(ij\right)-{\mathrm{tci}}_{\mathrm{NB}}\left(ij\right),\end{array}$$

where T_NB(ij) is the original temperature measured and $T_{\mathrm{NB}}^{\prime }\left(ij\right)$ is the corrected temperature (Fig. 7d).

It is important to remark that at this step it does not matter if the resulting image $T_{\mathrm{NB}}^{\prime }$ is correct in temperature because the temperature calibration of the NB is performed later (Sect. 4.3); so, for example, using the minimum value of the black-target image instead of the mean one, it would not change the final result, that is, the SO₂ content.

The blackbody measurements were usually performed at the start and after about 1 h of scene acquisition. This simple geometric correction does not always give optimal results; in some cases, the “ghost” effect is only partially removed, and some noise effects remain in the scene images. This is probably due to the different temperatures of the black target and the collected scene, as well as to any variations in the internal temperature of the instrument. Generally, nighttime measurements are less affected by these errors.

4.2 Matching MODTRAN with broadband clear-sky temperatures

In principle, the BB temperatures do not have to be corrected because of the factory calibration of the sensor (the digital counts-to-temperature conversion). Currently, uncooled microbolometer IR cameras with detection sensitivity between 8 and 14 µm are optimised for temperatures between 270 and 370 K, which are higher than the average clear-sky and most of the clouds' temperatures. The temperature accuracy guaranteed by the manufacturer is the larger of ±5 °C or ±5 % between 5–140 °C and ±10 % between 140 and 330 °C (https://www.thermal.com/uploads/1/0/1/3/101388544/mosaic_core_specification_sheet_2021v2-web.pdf, last access: 7 January 2025). The calibration of the BB channel was already confirmed in the laboratory between 270 and 300 K using a custom-built blackbody plate of high emissivity (ϵ > 0.98) that could be heated and cooled using a thermoelectric (Peltier) cooler that operates by applying a current to a semiconductor material (Prata et al., 2024). Conversely, at lower temperatures, the calibration is not optimal and the temperature of the chalcogenide lens (Fig. 1b), placed in front of the sensor, becomes predominant. Moreover, the lower limit detectable for this type of sensor is 230–240 K, while the clear-sky temperatures could easily reach lower values, especially at high altitudes, where generally the volcano measurements are carried out. The minimum temperature we measured was T_BB = 230.6 K at Sabancaya Volcano (Perù) at 5000 m a.s.l.

To make the quantitative retrieval effective, we used MODTRAN 5.3 to simulate the temperature measured by the TIR cameras. Even if a proper atmospheric pressure, temperature, and humidity (PTH) vertical profile is used (near in space and time to the VIRSO2 acquisitions), the scene collected by TIR cameras represents the real local atmospheric situation in a volcanic environment, with the presence of gases and particles (ash and water droplets) in an atmosphere which is not simple to simulate with a radiative transfer model (RTM). This, together with the low-temperature limits of BB sensor, very often causes the clear-sky T_BB to be higher than those temperatures simulated by the model. In Fig. 8a the green crosses represent the minimum temperatures computed row by row from T_BB for the scene of Fig. 7a. To avoid problems related to pixels affected by the plume and too close to the land (with higher sky temperatures due to heating of the ground), only the lines from 1 to 150 were considered. The blue line represents the simulated MODTRAN clear-sky temperatures, as a function of the elevation angle θ (Eq. 1). The atmospheric PTH vertical profile on 30 August 2024 at 06:00 UTC (37.5° N, 15.0° E) from the National Centers for Environmental Prediction (NCEP) dataset (Kalnay et al., 1996) was used in the model, and the top of the atmosphere (at 100 km) was considered non-emitting (T_s = 0 K). The clear-sky T_BB values are about 10–20 K higher than the simulated temperatures. To correct this gap, it is not possible to simply add or subtract a certain value (depending on θ) to all the pixels of the image. The plume pixels are characterised by higher temperatures than those of clear sky; depending on their opacity, the difference with MODTRAN can be smaller or even absent.

Figure 8 Etna-Montagnola on 30 August 2024, 05:22 UTC, clear-sky test case. (a) Comparison between MODTRAN BB simulations with T_s = 0 K (blue line) and T_s = 225.8 K (red line) and VIRSO2 T_BB (green crosses). (b) Comparison between MODTRAN BB (red line) and NB (black line) simulated temperatures with the corresponding VIRSO2 measurements. The blue and cyan crosses represent the NB temperatures before and after the calibration procedure respectively.

It is possible to modify T_BB below a certain temperature, making some reasonable hypotheses about the factory calibration and the characteristics of the chalcogenide lens or, alternatively, bring MODTRAN simulation data closer to clear-sky measurements, adding aerosol or water vapour to the simulated atmosphere, even though this is not always enough. A simple and valid method to match the MODTRAN simulation with the T_BB image is instead to consider the top of the atmosphere an emitting layer above the plume. This layer may represent either the sensor's inability to go below certain temperatures or the effect of high cirrus clouds (or both). The presence of the layer at temperature T_s does not imply a generic additive term which could lead to erroneous results, such as that due to temperatures higher than those of the air for extremely opaque plumes. Indeed, an increase in the columnar contents (ash or SO₂) of the plume progressively reduces the effect of the presence of the T_s layer. In this case, the global clear-sky radiance L_s as a function of θ that reaches the camera is

$$\begin{array}{}\text{(15)}& L_{\mathrm{s}}\left(\mathit{\theta }\right)={\mathit{\epsilon }}_{\mathrm{s}}\cdot B\left(T_{\mathrm{s}}\right)\cdot {\mathit{\tau }}_{\mathrm{0}}\left(\mathit{\theta }\right)+L_{\mathrm{0}}\left(\mathit{\theta }\right),\end{array}$$

where B is the Planck function, ε_s and T_s are the emissivity and the temperature of the sky at the top of the atmosphere, and τ₀ and L₀ are the atmospheric clear-sky transmittance and path radiance (obtained from MODTRAN). Assuming that there is at least 1 pixel inside the 320 × 240 array unaffected by the plume (clear sky), it is possible to compute the clear-sky radiance L_s by applying the Planck formula to the absolute minimum of T_BB(ij). Knowing the specific elevation angle of the minimum used, now the T_s value, needed for matching the simulated and the measured BB temperatures, can be easily computed from Eq. (15), assuming ε_s = 1 for the sky (blackbody). Finally, the overall MODTRAN look-up tables (all angles, all different plume columnar contents, both BB and NB channels; see the detailed description in Sect. 5) can be obtained from Eq. (15) using the computed temperature T_s and the specific overall transmittance and path radiance. In Fig. 8a the red line represents the MODTRAN clear-sky temperatures computed with T_s = 226 K, matched to the minimum T_BB. Note the border effect of BB data that, in this case, produces an unrealistic increase (of about 1–2 K) in the clear-sky temperature for angles greater than 40°.

Although it involves a longer calculation time, it is advisable to carry out this procedure (T_s estimation and look-up table (LUT) recomputation) for every VIRSO2 image to take into account possible changes to the local atmospheric situation and/or instrument's internal temperature variations. In Fig. 8b the comparison between BB and NB MODTRAN clear-sky simulations is shown, together with the corresponding VIRSO2 BB and NB measurements. The cyan cross symbols represent the NB temperatures after the final recalibration (see Sect. 4.3).

4.3 Narrowband calibration

As said in the introduction of Sect. 4, the presence of the 8.7 µm filter placed in front of the camera strongly affects the measurements, destroying any previous factory calibration of the broadband sensor. In Fig. 7a, T_BB min/max variation is more than 40 K, ranging from about 242 K (sky) to 287 K (land), while $T_{\mathrm{NB}}^{\prime }$ (the NB temperature after the “ghost” effect removal, Fig. 7d) ranges only from about 285 K (sky) to 290 K (land). The very limited $T_{\mathrm{NB}}^{\prime }$ variation (only about 5 K) makes the calibration particularly important and critical at the same time. The need for this second correction is also evident from the comparison between clear-sky $T_{\mathrm{NB}}^{\prime }$ and MODTRAN (blue crosses and black line in Fig. 8b).

In principle, the NB measurements can be written as

$$\begin{array}{}\text{(16)}& B\left(T_{\mathrm{NB}}^{\prime }\right)={\mathit{\tau }}_{\mathrm{f}}\cdot B\left(T_{\mathrm{NB}}^{\prime \prime }\right)+\left(\mathrm{1}-{\mathit{\tau }}_{\mathrm{f}}\right)\cdot B\left(T_{\mathrm{f}}\right),\end{array}$$

where τ_f and T_f are the transmittance and the temperature of the 8.7 µm filter and $T_{\mathrm{NB}}^{\prime \prime }$ is the temperature of the target, i.e. the expected recalibrated NB temperature. While τ_f is approximately known (about 0.8), the T_f is unknown and variable: the filter is outside of the VIRSO2 system and therefore exposed to air, wind, and sun, but at the same time, it is near to the camera, so it is affected by internal temperature variations too. Moreover, to use it for recalibration, T_f must be known with extremely good precision, given that usually it gives the greater contribution to Eq. (16), with T_f ≫ $T_{\mathrm{NB}}^{\prime \prime }$. For this reason, we decided to use a similar, but easier, numerical recalibration method based, in first approximation, on a simple linear relationship:

$$\begin{array}{}\text{(17)}& T_{\mathrm{NB}}^{\prime \prime }(i,j)=a\cdot T_{\mathrm{NB}}^{\prime }(i,j)+b.\end{array}$$

Given that there are two coefficients to be derived (a and b), a minimum of two different calibrated temperatures $T_{\mathrm{NB}}^{\prime \prime }$ are needed, preferably covering a wide range of temperatures, including those of the volcanic plume. A good choice for the lower one is the minimum sky temperature already used in the previous section, Sect. 4.2 (yellow box in Fig. 9a). For this point, we already know the expected value for the NB channel from the MODTRAN LUTs recomputed with T_s = 226 K. Generally, the simulated NB clear-sky temperature is lower than the BB temperature of about 5–10 K, depending on the atmospheric composition, the ground altitude, and the elevation angle of the camera (see the red and black lines in Fig. 8b).

Figure 9 Etna-Montagnola on 30 August 2024, 05:22 UTC, test case. (a) Calibrated and co-registered $T_{\mathrm{NB}}^{\prime \prime }$ temperature (the white and yellow boxes indicate the pixels considered in the ground–sky calibration). (b) Ground simulated δT_G values for a target placed at different distances (1–9 km) from the camera (θ = 10°). The target has a uniform emissivity (ε_G = 0.95) and a temperature range of 270–310 K.

On the other hand, a good candidate for the higher temperature value can be the average of a small ground area (white box in Fig. 9a). Land temperatures are indeed usually higher than those of the plume and, by choosing the area far enough from the crater, it will not be affected by SO₂. Specific MODTRAN simulations computed considering a “grey” target allow us to obtain the theoretical relationship between the BB and NB ground temperatures, which depends mainly on the ground temperature and the distance of the volcano from the camera, as well as on the atmospheric composition (Fig. 9b).

Assuming T_BB is correct, the calibrated $T_{\mathrm{NB}}^{\prime \prime }$ values for sky and ground are

$$\begin{array}{}\text{(18)}& {\displaystyle }{T}_{\mathrm{NB}}^{\prime \prime }\left(\text{sky}\right)=T_{\mathrm{BB}}\left(\text{sky}\right)-\mathit{\delta }{T}_{\mathrm{S}},\text{(19)}& {\displaystyle }{T}_{\mathrm{NB}}^{\prime \prime }\left(\text{ground}\right)=T_{\mathrm{BB}}\left(\text{ground}\right)-\mathit{\delta }{T}_{\mathrm{G}},\end{array}$$

where δT_S and δT_G are the sky and the ground temperature values to be subtracted from the BB to obtain the corrected NB temperature. The labels (sky) and (ground) represent the average temperatures of two small areas of the image related to clear sky and the ground respectively (Fig. 9a). As already reported in the introduction of Sect. 4, we compare BB and NB pixels in the same scene position as where the NB array must have been previously co-registered to BB. In the specific test case presented (Etna-Montagnola, 30 August 2024), with the VIRSO2 camera placed at 2600 m a.s.l. and at 3.5 km from the crater, δT_S = 6.2 K and δT_G = −0.1 K, from which the linear calibration coefficients of Eq. (17) become a = 9.76 and b = −2542. Figure 9a shows the resulting calibrated NB image, now ranging from 236 to 290 K.

4.4 “Zero image” and “temperature difference” computation

Once the recalibration and the other correction procedures have been applied, a temperature difference image is made, consisting of the difference between the NB-corrected image and a “zero value” corresponding to imagery with no plume (Eq. 20):

$$\begin{array}{}\text{(20)}& {\text{DT}}_{\mathrm{NB}}(i,j)=T_{\mathrm{NB}}^{\prime \prime }(i,j)-T_{\mathrm{0},\mathrm{NB}}^{\prime \prime }\left(i\right).\end{array}$$

A similar temperature difference image can be obtained from the BB camera too (Eq. 21):

$$\begin{array}{}\text{(21)}& {\text{DT}}_{\mathrm{BB}}(i,j)=T_{\mathrm{BB}}(i,j)-T_{\mathrm{0},\mathrm{BB}}\left(i\right).\end{array}$$

These images (Fig. 10), labelled DT_BB and DT_NB, are the main information content of the retrieval scheme; in particular, the DT_NB image will be used to obtain the SO₂ plume columnar contents by comparing it with the simulated MODTRAN values.

Figure 10 (a) BB and (b) NB temperature difference images (Etna-Montagnola, 30 August 2024, 05:22 UTC).

As already noted in Fig. 8a, some small border effects (vignetting) are present in the BB image (Fig. 10a), so we obtained some DT_BB values > 0 which are not related to the volcanic plume's presence. Obtaining the zero images is not so easy in practice as they should contain no plume and yet image the area of the sky close to where there is a plume so that a similar part of the atmosphere is sampled. Moreover, because the clear air temperatures are also not uniform and increase as the θ angle decreases (see Fig. 8), it is necessary to consider a “zero value” for each row of the image. It is possible to simply consider the smoothed row-by-row minimum values of the overall image or the row-by-row average values of a specific plume-free sub-region of the image. Alternatively, it is also possible to use the clear-sky temperatures simulated with MODTRAN after the “matching” procedure (red and black lines of Fig. 8b). These tend to be less noisy, but on the other hand they can produce some small biases in the DT values if the matching VIRSO2-MODTRAN is not accurate. The SO₂ retrievals are very sensitive to the variations in DT_NB, so an imprecise “zero image” computation can cause large retrieval errors.

5 SO₂ retrieval

Using MODTRAN, we can simulate the behaviour of BB and NB temperatures according to different SO₂ vertical columnar densities and elevation angles. This is done by inserting a horizontal layer with increasing SO₂ vertical column density (VCD) contents into the atmospheric vertical profile. This means that, for MODTRAN, the volcanic plume occupies an infinite (in all directions) horizontal layer of atmosphere. This configuration (“MODTRAN plume” in Fig. 11) would be correct only if the plume moved exactly in the direction of the camera (ω = 90°). If instead the wind direction is parallel to the camera's focal plane (ω = 0°), to obtain a correct configuration of the simulation, it is necessary to make a conversion from the “real plume” to the “MODTRAN plume” geometry (from vertical to horizontal). In Fig. 11, the scheme of the problem is shown. As θ increases, it is necessary to increase the height and the thickness of the MODTRAN layer where the plume is present to make sure the path length of the radiance is correct (see the segments CP₁, CP₂, P₁Q₁, and P₂Q₂ in Fig. 11, which are the same for both the configurations). In the MODTRAN configuration, the model does not account for upwelling and downwelling radiance from the lower and upper part of the plume, regarding a certain optical path. This is partially compensated by the presence of the plume on the left and on the right, regarding a certain optical path. In any case the multiple-scattering contribution (used in the MODTRAN simulations) is very low.

Figure 11 Example scheme of the real vs. MODTRAN plume configuration if the wind direction is parallel to the camera's focal plane. The VIRSO2 camera is at C, and the focal plane is perpendicular to the figure.

The following simple equations were used to compute the bottom and top altitude of the MODTRAN plume layer:

$$\begin{array}{}\text{(22)}& {\displaystyle }\text{BH}=\left(D_{\mathrm{0}}-{\displaystyle \frac{\text{HT}}{\mathrm{2}}}\right)\cdot \tan \left(\mathit{\theta }\right)+h_{\mathrm{0}},\text{(23)}& {\displaystyle }\text{VT}=\text{HT}\cdot \tan \left(\mathit{\theta }\right),\end{array}$$

where h₀ is ground altitude, BH is the bottom height of the plume (a.s.l.), HT and VT are the horizontal and vertical thickness of the plume, and D₀ is the horizontal distance to the target. If the wind direction is not parallel to the focal plane, that is, ω ≠ 0°, as a first approximation, D₀ can be simply replaced by the average $D_{\mathrm{0}}^{*}$ value (see Sect. 3.1). HT is an input parameter which is usually unknown, but it does not significantly affect the results (see Sect. 7).

Figure 12 shows the simulated DT_BB and DT_NB vs. SO₂ VCD, varying the viewing angles θ from 10 to 50°, for the Etna-Montagnola test case (30 August 2024, 06:00 UTC, NCEP PTH profile; h₀ = 2600 m, $D_{\mathrm{0}}^{*}$ = 3 km, and HT = 200 m). As mentioned in Sect. 4.2, the MODTRAN LUTs were recomputed with a more appropriate sky temperature, in this case, T_s = 226 K. It is clear a strong dependence of DT on elevation angles. In case of presence of SO₂ only (no particles in the plume), in the BB image the plume should be poorly visible as simulated DT_BB ranges only between 0 and 5 K. This is not the case for Fig. 10a, a clear sign that in the depicted emission there are also some particles (ash and/or water droplets). The NB channel is instead much more sensitive to SO₂, with MODTRAN DT_NB ranging between 0 and 32 K. Knowing θ(ij) for each pixel of the image, by means of a simple bilinear interpolation, the final image of SO₂ column content can then be computed from the DT_NB image. Figure 13a shows the SO₂ VCD obtained as described, which attributes all the NB signal (DT_NB) to SO₂ only. Note that to mitigate the inaccurate “ghost” image removal (Fig. 10b), only the pixels with DT_BB > 2 K were considered. In the upper part of the image, some edge effects still remain. These pixels must be excluded in the computation of SO₂ mass and flux (Sect. 6).

Figure 12 BB (a) and NB (b) MODTRAN DT vs. SO₂ VCD performed for the Etna-Montagnola test case (30 August 2024).

Figure 13 Etna-Montagnola on 30 August 2024. (a) A 05:22 UTC SO₂ VCD image (g m⁻²) obtained from the NB channel, considering SO₂ only (not corrected for the particles). (b) A picture of the volcanic plume at 05:00 UTC.

Here the presence of particles in the plume (ash and/or water droplets) was not considered, while such particles strongly affect both BB and NB channels. Figure 13b shows a close-in-time picture of the volcanic plume in which the presence of condensed water vapour particles seems highly probable. The simplification to attribute all the NB signal (DT_NB) to SO₂ only produces an overestimation of the results obtained. Removing the contribution of the particles from DT_NB will be an important improvement of the algorithm in the future.

6 SO₂ mass and flux computation

The comparison with MODTRAN gives, for each pixel of the image, the SO₂ VCD content. The computation of SO₂ mass (M_s) and flux (F_s) requires converting the “MODTRAN” vertical content to the “real” slant column density (SCD, Eq. 24 and Fig. 11). Now the SO₂ total emitted mass can be simply computed multiplying the SCD for the area of the pixels (Eq. 25; Tamburello et al., 2012).

$$\begin{array}{}\text{(24)}& {\displaystyle }\text{SCD}=\text{VCD}/\sin \left(\mathit{\theta }\right)\text{(25)}& {\displaystyle }{M}_{\mathrm{s}}={\sum }_{i,j}a(i,j)\cdot \text{SCD}(i,j)\end{array}$$

In this way, it is possible, for example, to compare VIRSO2 results with other ground-based sensors such as UV cameras, paying attention to ensuring that the same FOV area is considered. Instead, to validate the VIRSO2 results with space-based data (satellite), because of the very different points of observation and spatial/temporal resolutions, the SO₂ flux is a better quantity to use. Considering a transect perpendicular to the direction of the plume (vertical, horizontal, or oblique), the SO₂ flux (function of time) is computed:

$$\begin{array}{}\text{(26)}& F_{\mathrm{s}}\left(t\right)={\sum }_{i=\mathrm{1}}^N{v}_i\left(t\right)\cdot {\text{SCD}}_i\left(t\right)\cdot l_i,\end{array}$$

where i (from 1 to N) denotes the pixels contained in the transect, v is the plume velocity, and l is the pixel width (orthogonal to the flux direction). The plume velocity can be computed pixel by pixel and time step by time step using an optical flow algorithm or simply obtained from the focal plane component of the wind speed (ws) and considered the same for all the plume pixels (Eq. 27):

$$\begin{array}{}\text{(27)}& v=\text{ws}\cdot \cos \left(\mathit{\omega }\right).\end{array}$$

In this case (Etna-Montagnola, 30 August 2024) the volcanic emission is quite low, and the optical flow does not give reliable results. So, the wind speed was set to 2.1 m s⁻¹, obtained as the average speed at 3.5 km a.s.l. at Etna (37.5° N, 15.0° E) on 30 August 2024 at 05:00–07:00 UTC from the ERA5 hourly dataset. In the same way, the average wind direction was set to 38° north, so ω ≈ 30° and the plume velocity becomes v = 1.8 m s⁻¹. Finally, it is recommended to choose a transect not too close to the crater because in the RTM simulations the plume is assumed to be in thermal equilibrium with the surrounding air temperature. For this reason, the VIRSO2 flux is computed as the average of the fluxes obtained considering the 75th, 100th, and 125th columns as vertical transects. The first 40 lines were discarded due to some noise/residual error (ghost effect) and because the plume generally stayed in the lower part of the images for this test case (see Fig. 13a and BB temperature videos in the Supplement).

The TROPOspheric Monitoring Instrument (TROPOMI) is a spectrometer on board the Sentinel-5 Precursor (S5P) polar orbit satellite that covers a spectral range from ultraviolet (UV) to shortwave infrared (SWIR), with a spatial resolution at the nadir of 5.5 × 3.5 km² and a revisit time of about 1 d (Veefkind et al., 2012). Particularly significant for TROPOMI UV bands is the ability to detect SO₂ columnar abundance of about 0.02 g m⁻² (0.7 DU), i.e. sensitivity about 30 times better than the sensitivity of the multispectral satellite sensors (Corradini et al., 2021). Here the near-real-time (NRT) Level 2 SO₂ product collected on 30 August 2024 at 11:33 UTC is considered. The SO₂ VCD image (Fig. 14a) is obtained by interpolating at 3.5 km a.s.l. the two “sulfur dioxide total vertical column” products at 0.5 km a.g.l. and 7 km a.s.l. Finally, the SO₂ flux is computed by applying the traverse method, and the wind speed is the same as that considered for VIRSO2 (ws = 2.1 m s⁻¹). Differently from the camera, all the transects perpendicular to the volcanic cloud axis are considered (Pugnaghi et al., 2006; Merucci et al., 2011; Theys et al., 2013; Corradini et al., 2021). Depending on wind speed and plume extension, this makes it possible to obtain an SO₂ flux trend for several hours before the image acquisition (shown partially in Fig. 14b).

Figure 14 Etna-Montagnola on 30 August 2024. (a) The S5P-TROPOMI SO₂ VCD map (11:33 UTC) from which the SO₂ flux has been computed. (b) SO₂ flux (in metric tons per day) obtained from TROPOMI (orange line and symbols) and from VIRSO2 (cyan symbols and blue lines). The solid blue line is a moving average over 1 min, while the shading represents an error of ±40 % in the retrieval.

Figure 14b shows the SO₂ fluxes obtained from the TIR camera and TROPOMI, with the first one greater than about 300–800 t d⁻¹ (metric tons per day). The mean values are 481 and 996 t d⁻¹ for TROPOMI and VIRSO2 respectively. The VIRSO2 retrieval clearly shows the fluctuations in emission, which are not visible from satellite. Considering a ±40 % error for both VIRSO2 (see Sect. 7) and TROPOMI (±35 % retrieval error (Corradini et al., 2021) ±20 % for wind speed, quadrature sum), the data mainly fall within the uncertainty of the procedure. In this case, as said in Sect. 5, the overestimation of VIRSO2 is probably due to the presence of water vapour particles condensed in the volcanic plume.

7 Sensitivity study

Due to the various steps and inputs of the procedure, there are several possible sources of error, not all of which are easily quantifiable. For example, the ghost image removal, the NB channel recalibration, and the “zero” image computation can introduce errors that cannot be assessed a priori and that can have an impact on the SO₂ VCD retrieval. Other sources of errors, related to the uncertainties of the input parameters, are more easily computable.

Table 2 Etna-Montagnola on 30 August 2024 (2600 m a.s.l.). SO₂ flux errors related to the input parameters of the procedure.

The sensitivity study was done by evaluating the total mass that erupted on 30 August 2024 between 05:30 and 07:18 UTC and was obtained by integration of the SO₂ flux. For the test case considered, the reference configuration, its variability and the related errors are shown in Table 2. The impact on SO₂ of any errors in the MODTRAN simulations was evaluated through the first three input parameters (HT, D₀, and PTH). As Table 2 shows, while the geometrical configuration (HT and D₀) uncertainties do not lead to significant retrieval errors, the use of an appropriate atmospheric vertical profile (PTH) is important. In this case, the August and September 1981–2010 monthly mean profiles around Etna (32.5–42.5° N, 10.0–20.0° E) were considered instead of the specific date (30 August 2024, 06:00 UTC), obtaining an error of −7.6 % and 14.2 % respectively. The central elevation angle (α) and the wind direction mainly affect the size of the pixels of the image, which is important for the computation of both SO₂ total mass and flux (Eqs. 25–26). α is a very sensitive parameter because it affects the SCD computation too (Eq. 24), but it can usually be obtained with good precision using an inclinometer attached to the VIRSO2 camera or from the altitude of a known reference point in the image, for example the top of the volcano (Snee et al., 2023). The same accuracy cannot be achieved for the wind direction. This latter parameter can be obtained from reanalysis data (as in this case) or directly from in situ measurements, but the local effects related to orography can produce significant variations. Finally, the plume speed is directly proportional to the flux (see Eq. 26), so its uncertainty proportionally affects the flux.

In summary, an overall uncertainty of about ±40 % can be calculated as the sum in quadrature of the maximum (as an absolute value) individual uncertainties.

8 Conclusions and future work

In this paper, we present a new TIR ground-based system for the detection of volcanic plumes and the retrieval of geometry and SO₂ emissions. The main advantages of this instrument are its high portability, its relatively low cost, and the possibility of making measurements during both the day and the night. The drawbacks are mainly related to the calibration: the presence of the filter centred at 8.7 µm outside of the system produces some critical effects that need to be corrected through a post-processing procedure. The overall retrieval procedure consists of three different steps: first the instrument view geometry is computed, then the calibration is carried out, and finally the SO₂ is retrieved. Every step of the process introduces errors: here an overall error of ±40 % was estimated for the SO₂ flux retrievals by considering the uncertainties in the instrument viewing geometry and the atmospheric profiles of pressure, temperature, and humidity.

As test cases, some VIRSO2 measurements were performed at Mt. Etna in March–April 2021 and August 2024. The measurements collected in 2021 show the effectiveness of the plume detection and geometry retrievals and the importance of the wind correction, while the 2024 SO₂ flux retrievals show results comparable with those obtained from space using the S5P-TROPOMI UV sensor despite the huge difference between the two instruments (spectral range, spatial resolution, point of view). In particular, the higher values of SO₂ flux obtained from the camera compared with TROPOMI could be due to the presence of particles (ash or water droplets) in the volcanic plume that cause an overestimation of the VIRSO2 retrieval. The aim of the article is mainly methodological, that is, to explain in detail the algorithm developed and used for the SO₂ estimation. The comparison and validation with other instruments, both satellite and ground-based, will be the subject of other papers (in preparation) in which the data collected at Etna, Stromboli, Sabancaya (Peru), and Popocatepetl (Mexico) volcanoes will be used.

Several improvements can be included in the procedure in the future. A low-temperature blackbody instrument, composed of a uniform emissive head, a refrigerating unit, and an electronic temperature controller, would be very useful to understand the behaviour of the BB sensor below 270 K, which is outside of the factory calibration range of the camera. The ghost image removal procedure presented in this paper is very simple but sometimes not completely effective due to the temperature difference between the black target and the real detected scene. The incomplete elimination of the “ghost” effects can produce errors in the following recalibration of the NB 8.7 µm channel, which is the crucial part of the algorithm, so an improvement of this part would be desirable. The addition of temperature sensors near the filter and inside the camera body would help and be quite easy to implement. Finally, the presence of particles in the volcanic plume (ash, water droplets, ice, etc.) strongly affects the NB 8.7 µm measurements. The correction for particles' contribution will represent an important improvement of the algorithm.

The SO₂ flux is a fundamental parameter for the study of the volcano's activity and for comparison with other instruments, from both ground and space. To compute it, the velocity of the plume is a necessary value to know; in this sense, an optical flow algorithm could give important information regarding the different speeds of the various parts of the plume, in both time and space.