# Plume Height Time-Series Retrieval Using Shadow in Single Spatial Resolution Satellite Images

Sophie Pailot-Bonnétat, Andrew J L Harris, Sonia Calvari, Marcello de Michele, Lucia Gurioli

► To cite this version:

Sophie Pailot-Bonnétat, Andrew J L Harris, Sonia Calvari, Marcello de Michele, Lucia Gurioli. Plume Height Time-Series Retrieval Using Shadow in Single Spatial Resolution Satellite Images. Remote Sensing, 2020, 12 (23), pp.3951. 10.3390/rs12233951. hal-03053523

HAL Id: hal-03053523

https://uca.hal.science/hal-03053523

Submitted on 11 Dec 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Article

# Plume Height Time-Series Retrieval Using Shadow in Single Spatial Resolution Satellite Images

Sophie Pailot-Bonnétat <fnref n="1" />*, Andrew J. L. Harris <fnref n="1" />, Sonia Calvari <fnref n="2" />, Marcello De Michele <fnref n="3" /> and Lucia Gurioli <fnref n="1" />

Received: 13 October 2020; Accepted: 30 November 2020; Published: 3 December 2020

<figure><img src="image_3.png" /></figure>

check for
updates

**Abstract:** Volcanic plume height is a key parameter in retrieving plume ascent and dispersal dynamics, as well as eruption intensity; all of which are crucial for assessing hazards to aircraft operations. One way to retrieve cloud height is the shadow technique. This uses shadows cast on the ground and the sun geometry to calculate cloud height. This technique has, however, not been frequently used, especially not with high-spatial resolution (30 m pixel) satellite data. On 26 October 2013, Mt Etna (Sicily, Italy) produced a lava fountain feeding an ash plume that drifted SW and through the approach routes to Catania international airport. We compared the proximal plume height time-series obtained from fixed monitoring cameras with data retrieved from a Landsat-8 Operational Land Imager image, with results being in good agreement. The application of the shadow technique to a single high-spatial resolution image allowed us to fully document the ascent and dispersion history of the plume-cloud system. We managed to do this over a distance of 60 km and a time period of 50 min, with a precision of a few seconds and vertical error on plume altitude of ±200 m. We converted height with distance to height with time using the plume dispersion velocity, defining a bent-over plume that settled to a neutral buoyancy level with distance. Potentially, the shadow technique defined here allows downwind plume height profiles and mass discharge rate time series to be built over distances of up to 260 km and periods of 24 h, depending on vent location in the image, wind speed, and direction.

**Keywords:** Etna volcano; lava fountain; paroxysmal explosive eruptions; ash plume height; Landsat 8 satellite images; mass discharge rate time-series

## 1. Introduction

Volcanic ash clouds injected into the atmosphere represent a threat to population health and aircraft operations [1]. This is a recurrent problem, and despite the improved monitoring of volcanic clouds, aircraft encounters and incidents occur regularly (e.g., [1–5]). Damage to aircraft, traffic rerouting and airport traffic disruption can cause millions of dollars in losses for airline companies and any related businesses [6]. Even explosive eruptions at volcanoes located in unpopulated areas represent a threat because the associated ash clouds can be transported thousands of kilometers away from their sources [7]. In volcanic hazard assessment, the retrieval of cloud height is a key parameter because it can be linked to dispersion [8,9], and is used to derive other, higher-level parameters that allow better classification of the eruption and understanding of the associated eruption dynamics, such as

<footnote n="1">¹ Laboratoire Magmas et Volcans, Université Clermont Auvergne, CNRS, IRD, OPGC, 63000 Clermont-Ferrand, France; andrew.harris@uca.fr (A.J.L.H.); lucia.gurioli@uca.fr (L.G.)</footnote>

<footnote n="2">² Osservatorio Etneo—ezione di Catania, Istituto Nazionale di Geofisica e Vulcanologia, Piazza Roma 2, 95125 Catania, Italy; sonia.calvari@ingv.it</footnote>

<footnote n="3">³ BRGM, Bureau de Recherches Géologiques et Minières, 3 Avenue Claude Guillemin, F-45060 Orléans, France; m.demichele@brgm.fr</footnote>

<footnote>* Correspondence: sophie.pailot_bonnetat@uca.fr; Tel.: +33-(0)631358090</footnote>
the Volcanic Explosive Index (VEI) [10] and mass discharge rate (MDR, i.e., mass flux per unit time in kg s⁻¹) [8].

The synoptic view offered by satellites is especially useful in observing large-scale phenomena such as volcanic clouds [11,12]. Other advantages of satellite data include their free or cheap availability, rapid (near-real-time) dissemination, synoptic view and global coverage [7]. Moreover, many volcanoes that experience low-frequency but high-magnitude and -intensity—Vulcanian through Plinian—eruptions [13,14], do not possess ground-based monitoring networks so that sometimes satellite data are the only available data to detect, track and study an eruption [11,12,15]. Satellites have thus been used for an increasing number of studies focused on detecting volcanic clouds and tracking their dispersal since the 1970s. As part of this, several methods that use satellite sensor-based data have been developed to retrieve cloud heights (as reviewed in Supplementary Materials Table S1). One of them is the height-from-shadow technique, whereby the length of the shadow cast by the cloud onto the ground is measured, and a “simple” trigonometric relationship with the sun geometry delivers the cloud height [16]. However, the shadow technique has not frequently been used for volcanological applications, especially using high-spatial resolution (30 m pixel) data, such as LANDSAT’s MSS, TM, ETM+ and OLI, Terra’s ASTER and ESA’s Sentinel (cf. Supplementary Materials Table S2). The height-from-shadow method does have some issues, such as the assumptions that need to be made regarding cloud shape or plume topography [17]. It is also only applicable during daytime when shadows are present. It is, however, unclear why such high-spatial resolution data have not been used to more frequently retrieve heights for volcanic clouds. This is especially so when considering problems associated with the commonly used height-from-temperature or brightness temperature method, which suffers from problems such as plume transparency or ice formation [7,18–21], and which are more typically applied to low-spatial resolution (1–4 km pixel) data from sensors such as AVHRR, MODIS, GOES, and SEVIRI. The shadow technique also has the advantage of not requiring stereo-imaging or temperature measurement, nor does it require assumptions regarding cloud thermal and/or emissive and transmissive properties. In this study, we thus develop and validate a method that allows the shadow technique to be used to (i) retrieve cloud height in high-spatial resolution data, and (ii) convert the downwind plume heights to at-vent MDR time series.

A distinction needs to be made between volcanic plumes, which are connected to the source (rooted to the vent), and volcanic clouds, which have become detached from the source [22]. The former is characterized by a dominant vertical component, the latter by a dominant horizontal component. Eruptive columns can be divided from bottom to top into three regions. The first and lowermost is the gas-thrust region immediately above the vent where the plume rises as a momentum-driven jet [8,23,24]. The second is the convective region, where the density of the plume becomes less than that of the surrounding atmosphere, atmospheric mixing occurs and ascent by buoyancy takes over [8,23]. The third is the umbrella region, where the densities of the plume and the atmosphere become equal and the cloud spreads laterally, defining the neutral buoyancy level (NBL) [23,25]. Initially, though, the maximum plume height will be higher than the NBL because of the inertia of the particles arriving at the top of the column to lead to overshooting [23]. Thus, while the dispersal of volcanic clouds depends on their interaction with atmospheric motions (mostly wind speed and direction), plume dynamics will also depend on atmospheric stratification through the entire height of their ascent, variation in wind speed with height, as well as “source” conditions, notably MDR and air entrainment rates [8,22]. Depending on the eruption intensity (MDR) and wind strength (speed), we can find a continuum of plume behaviors. The end members are umbrella clouds, which are strong plumes that rise above the vent and spread horizontally around the neutral buoyancy height, and bent-over or weak plumes [26] that are significantly distorted downwind and may not reach the NBL before being sheared-off downwind [8].

For this study, a cloud-free Landsat-8 Operational Land Imager (OLI) image of the 26 October 2013 lava fountaining event at Mount Etna (Sicily, Italy) was selected. This target has the advantage of possessing excellent ground-truth data, as are available from the continuously monitoring camera
network operated by the Istituto Nazionale di Geofisica e Vulcanologia (INGV—Osservatorio Etneo, Catania, Italy). This consists of a network of thermal and visible cameras (Figure 1), which have been previously applied to plume tracking [27] and which are routinely used for monitoring. Their specifications are detailed in Ref. [28], where the GPS-stamped image frames allow for precisely timed event detection. In processing the Landsat 8 data, to produce plume or cloud height profiles versus distance from the vent, several parameters such as plume topography, shape, ground topography and satellite viewing geometry need to be taken into account [cf. 17]. This is done here through developing a methodology to extract height with distance using the cloud shadow. We convert these profiles into times series using cloud dispersal velocities obtained from the INGV camera network. This allows us to plot plume height over the vent for one hour prior the satellite image acquisition and to derive an MDR time series. This permits us to (i) describe the effects of settling to the NBL, (ii) assess the effect of wind on the bending of a fire fountain, and (iii) derive physical and empirical relations for the ascent and dispersal of this fountain-fed plume-cloud system.

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Ortho-rectified Landsat 8 OLI image of Mt. Etna acquired on 26 October 2013 at 09:37 UTC (courtesy of National Aeronautics and Space Administration (NASA)/United States Geological Survey (USGS)). Flight routes and altitudes for air traffic approaching Catania Fontarossa airport are taken from www.flightradar24.com. The two cameras from the INGV—Osservatorio Etneo camera network used in this paper are marked. A zoom on each cloud is given as inset. SV = Simeto Valley; CTA = Catania Fontanarossa airport.</figcaption></figure>
## 2. Case Study: OLI Imagery of the 25–26 October 2013 Fountaining at Etna

Mount Etna is characterized by Strombolian and lava fountaining activity, and frequent effusive eruptions (e.g., [29,30]). The climaxes of lava fountaining at Etna, which gave rise to ash plumes and clouds, are termed “paroxysmal episodes” [28]. The fall out of lapilli and ash represents a frequent hazard to the population residing on the flanks of Etna [31], which is around one million [32]. The threat to Catania airport is also a persistent hazard [5] where, for example, between 2011 and 2015 alone there were 49 lava fountain events [28], and 150 between 1990 and 2015 [33]. Fountaining and Strombolian, as well as Vulcanian, activity can be sourced from flank vents feeding effusive activity, as well as from the summit craters. Etna hosts five active summit craters: Bocca Nuova (BN), Voragine (VOR), North-East Crater (NEC), South-East Crater (SEC) and the New South-East Crater (NSEC) [28,34]. Of these, SEC has become the most common locus of lava fountain activity since its formation in the 1970s (e.g., [28,35]). The NSEC developed from a pit on the east flank of the SEC that formed in 2011 and fed fountaining episodes between 2011 and 2015 [29,30]. It is one of these events, that of 25–26 October 2013, that we focus on here.

### 2.1. *The 25–26 October 2013 Eruption*

The 25–26 October 2013 eruption occurred after a six-month-long pause in activity [36]. It was marked by activity from several craters at the same time, with Strombolian explosions starting on 25 October at the NSEC. Its intensity increased during the night, and a lava fountain began around ~01:30 (hh:mm; all times are Coordinated Universal Time, UTC) on 26 October [36] feeding a light-toned ash plume. Lava flowed out of the NSEC at ~03:15 [36] and the NEC also started to emit a dark-gray ash plume at ~06:20 that dispersed at a lower altitude than the NSEC plume. At the same time, occasional explosions occurred at BN. The lava fountain at NSEC lasted until around 10:00, reaching an average height of 430 m [28]. Its plume consisted initially mostly of gas, but the amount of pyroclastics increased with the increase in intensity, causing tephra fallout. The plume was, though, rapidly bent-over toward the west southwest by the wind [37].

A satellite image from Landsat 8's OLI sensor was acquired at 09:37 on 26 October around the end of the lava fountain episode (Figure 1). The two clouds from the NSEC and the NEC are distinguishable from the tonal difference and are visibly being blown toward the WSW by the wind. The light-toned cloud from the NSEC seems to be separated into two parts, the first extending 8 km from the vent; then a second, after a gap of 22 km, being at least 30 km long (it running off of the western edge of the image). The two parts will be referred to, hereafter, as cloud 1 for the near vent part and cloud 2 for the distal part. Cloud heights for this event are available for validation purposes from the INGV camera network (Figure 1), and have been derived from this same image using an independent, parallax-based method, i.e., the Plume Elevation Model (PEM) of De Michele et al. (2019) [38], thus also allowing inter-method comparison.

### 2.2. *Data*

Since the method developed here needs a clear view of the shadow cast by a plume, or cloud, on the ground, only satellite images taken during daylight can be considered. The meteorological cloud cover needs to be as little as possible to ensure the best visibility, even though shadows cast over meteorological clouds can be used, provided that the meteorological cloud top height is known or can be calculated at the same time. The volcanic cloud must also be sufficiently opaque so as to provide a clearly defined shadow. However, even if the cloud is not dense or is partially transparent, the method can be applied as long as a shadow is visible. This is an advantage over the brightness temperature method or the PEM procedure, where transparency or lack of optical thickness can cause plume height underestimations [7,38].

The LANDSAT 8 OLI image was downloaded directly from the Earthdata Search engine (https://search.earthdata.nasa.gov/search), which provides access to NASA's Earth Observing System Data
and Information System (EOSDIS) services. The product available is a level 1 product, which is terrain-corrected. The pre-processing of such level 1 products includes geo-referencing (alignment of imagery to its correct geographic location) and ortho-rectification (correction for the effects of relief and view direction on pixel location) to ensure the exact positioning of the image [39]. For a cloud, though, this is not ideal, as geo-referencing can introduce distortions when a cloud is projected onto the surface [40]. LANDSAT 8's OLI and Thermal Infrared Sensor (TIRS) provides 11 different bands of data at a spatial resolution of 15–100 m (Table 1).

<table><caption>**Table 1.** LANDSAT 8's OLI and TIRS instrument characteristics.</caption><thead><tr><th>Bands</th><th>Wavelength (μm)</th><th>Spatial Resolution (m)</th><th>Measure</th></tr></thead><tbody><tr><td>Band 1—Coastal aerosol</td><td>0.43–0.45</td><td>30</td><td>TOA/SR</td></tr><tr><td>Band 2—Blue</td><td>0.45–0.51</td><td>30</td><td>TOA/SR</td></tr><tr><td>Band 3—Green</td><td>0.53–0.59</td><td>30</td><td>TOA/SR</td></tr><tr><td>Band 4—Red</td><td>0.64–0.67</td><td>30</td><td>TOA/SR</td></tr><tr><td>Band 5—NIR</td><td>0.85–0.88</td><td>30</td><td>TOA/SR</td></tr><tr><td>Band 6—SWIR 1</td><td>1.57–1.65</td><td>30</td><td>TOA/SR</td></tr><tr><td>Band 7—SWIR 2</td><td>2.11–2.29</td><td>30</td><td>TOA/SR</td></tr><tr><td>Band 8—Panchromatic</td><td>0.50–0.68</td><td>15</td><td>TOA</td></tr><tr><td>Band 9—Cirrus</td><td>1.36–1.38</td><td>30</td><td>TOA</td></tr><tr><td>Band 10—Thermal Infrared (TIRS) 1</td><td>10.6–11.19</td><td>100</td><td>Temperature</td></tr><tr><td>Band 11—Thermal Infrared (TIRS) 2</td><td>11.50–12.51</td><td>100</td><td>Temperature</td></tr></tbody></table>

NIR = near infrared; SR = surface reflectance; SWIR = shortwave infrared; TIRS = thermal infrared; TOA = top of atmosphere reflectance.

# 3. Methodology

## 3.1. *Height-from-Shadow Time Series from a Single High-Spatial Resolution Image*

The height-from-shadow technique is based on the fact that the plume and cloud casts a shadow on the ground. If the length of the shadow and the angle between the ground and the sunray direction is known, then basic geometry delivers an approximation of the plume height [17,19,41–46]. We here lay out a six-step methodology that begins with the definition of the fundamental Sun–Earth astronomical relationships that underpin the methodology, and ends with the generation of a cloud height time series.

### 3.1.1. Sun–Cloud Geometric Relations

The Sun–Earth astronomical relationships from Iqbal (1983) [47] were used to calculate the sun elevation and sun azimuth at the time of image acquisition, as converted to true solar time by taking into account the latitude and longitude of the summit of Etna. To facilitate the length measurements, the North-oriented image was rotated clockwise by the sun azimuth angle (Ψ) using cubic convolution. As a result, the sunrays are orientated in a vertical direction, so that shadow length along a correct orientation can be made along easy-to-define vertical lines. (Figure 2) A second line, called the intersection line, was composed of as many segments as there were direction changes in the cloud, and was drawn following the center of clouds 1 and 2. This line represents the position of the cloud height profile that will be created. The image was contrast-enhanced to maximize the intensity of the shadow darkness and the contrast with the background. As shown in Figure 2, a shadow length measurement (SLM) was carried out for each peak and trough in the enhanced shadow outline (i.e., at each point of variation, the shadow being the same between each turn-around point). Each selected point on the shadow outline was then linked to its matching point on the intersection line, following a vertical line, which now (following the rotation) represents the sun azimuth direction (Figure 2). The measurements were made in pixels and converted to distance (in m) using the pixel size calculated to take into account off-nadir and Earth curvature distortion following Harris (2013) [48].
Finally, each shadow length measurement was linked to a distance from the vent using the geographical coordinate of each pixel and taking into account Earth curvature.

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Schematic satellite view of the rotated plume/shadow system and parameter definition.</figcaption></figure>

### 3.1.2. Viewing and Shadow Geometry

Whether the entire length of the shadow is visible or not depends on the viewing geometry of the satellite with respect to the sun geometry (Figure 3a). Taking into account the distance between the Sun and the Earth relative to the size of the object that casts a shadow on the Earth's surface, the sunrays are assumed to be parallel lines. If the entire length of the shadow is visible, then a simple geometry is valid. However, if there is no area on the ground illuminated by the sun between the cloud and its shadow, then the distance needed for the calculation is missing a part of the shadow which is hidden behind the cloud (Figure 3b). In the case of a hidden shadow on a horizontal surface, the following correction can be applied:

$$ \text{Plume altitude (hidden shadow case)} = D' \tan \alpha + \frac{(D' \tan \alpha)^2}{D' \tan \sigma - D' \tan \alpha} \quad (1) $$

where $\alpha$ is the sun elevation and $\sigma$ the satellite scan angle. In our case, there was no hidden part, so Equation (1) did not need to be applied.

Four cases need to be considered depending on whether shadows are also apparent on the cloud top or not. In the first case, there are no shadows on the plume top, meaning that the plume top is relatively flat (Figure 4a). In the second case, there are shadows on the plume top, meaning that two different heights can be measured: the cloud hedge height, and the variation in height between the cloud edge and structures within the cloud (Figure 4b). In the third case, the shadow extends to the edge of the cloud, meaning only the cloud top height can be measured (Figure 4c). The final case is one intermediate between cases i and iii (Figure 4d). For the imaged case, the high-spatial resolution allowed us to determine visually that cloud 1 belongs to cases ii, iii and iv, and cloud 2 belongs to case i. As a result, our height measurements use the length of the cloud shadow on the ground (D) to obtain the height of the cloud edge (H), and the lengths of shadows cast on the cloud top (d) to obtain
variations in the height (h) across the cloud depending on the case identified (Figure 4). The case that needed to be applied to any point in the cloud was determined by visual inspection of the imagery.

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Satellite and Sun viewing geometry responsible for a "missing" shadow portion. (a) General geometry. (b) Detailed geometry for the case where a portion of the shadow is hidden.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Plume/cloud height measurement at nadir case depending on the shadow distribution. (a) Case i: flat plume. (b) Case ii: plume with topography where measurement of H and h is possible. (c) Case iii: plume with topography where measurement of H and h is possible. (d) Case iv: intermediate case: only the measurement of H+h is possible.</figcaption></figure>

### 3.1.3. Height Calculations and Corrections

The simple trigonometric relationship between sun elevation, shadow length and cloud height uses a distance projected onto a horizontal plane. However, the shadow length visible on the satellite image is modified by the slope of the terrain onto which it is projected. Two two-dimensional cases were considered. For case 1, the shadow is cast onto a surface tilted in the same direction as the sunrays (Figure 5a) and in the opposite direction in case 2 (Figure 5b). The plume altitude for each case can be calculated following:

$$ \text{Plume altitude (case 1)} = \tan \alpha \left( \frac{H1}{\tan(\sin^{-1}(\frac{H1}{SLM}))} \right) + AI \quad (2) $$

$$ \text{Plume altitude (case 2)} = \tan \alpha \left( \frac{H1}{\tan(\sin^{-1}(\frac{H1}{SLM}))} + \frac{H1}{\tan \alpha} \right) + AI \quad (3) $$

where H1 is the elevation difference between the cloud shadow edge and its matching point in the intersection line, SLM is the shadow length measurement, and AI is the altitude of the point on the intersection line. This last value is added to the above-surface height measurement to obtain the plume altitude in meters above sea level (m a.s.l.), rather than meters above the ground surface (Figure 5b). The sign of H1 determines the case: if H1 is negative, then the geometry can be described by case 1; if positive, by case 2. The ground surface altitudes were obtained using Google Earth Pro with a precision of ±5 m.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Effect of topography on shadow length measurement. (a) Case 1: surface tilted in the same direction as the sunray direction. (b) Case 2: surface tilted in the opposite direction to the sun ray direction.</figcaption></figure>

Next, we need to correct for the effect of the ortho-rectification, where the raw data are projected onto a digital elevation model. While this insures that all ground points are located correctly and image distortion is removed [39], it also projects airborne objects (i.e., clouds) from their position above the ground onto the topography so that they consequently become distorted. As a result, the shadow edge line is at its correct location on level 1 data, but the cloud edges and cloud-top features are not (Figure 6a). The cloud appears, instead, as a flat feature that is draped across the topography. This affects the shadow length measurements. We can obtain the shadow length measurement in geometrically corrected level 1 (SLM1) data, but we need the shadow length measurement in reality (SLM<sub>r</sub>), i.e., corrected for geometric distortion, if we are to correctly calculate the actual cloud height, h<sub>p</sub> (Figure 6b). In the case of a projection onto a plane in the direction of the sunrays, the problem can be resolved in two dimensions (Figure 6b). In Figure 6b, H1 is the altitude difference between plume edge and the shadow edge as projected onto the topography in the ortho-rectified image. H2 is the altitude difference between the cloud edge as located on the ortho-rectified image and the intersection with the sunrays above this point on the ground. This, effectively, projects the position of the cloud edge from its flattened position on the ground vertically upwards onto its actual above-ground position. As a result, we define parameter x, which is the altitude difference between the cloud edge as placed on the level 1 image and the actual location of the cloud edge in the air (Figure 6b). We thus first need to calculate SLM<sub>r</sub> (Equation (4)), which then allows us to estimate h<sub>p</sub> using Equation (5):

$$SLM_r = \frac{SLM_1 * [H1 + H2 + (\tan \alpha * \frac{H2*SLM_1}{\tan \sigma - H1 - H2})]}{H1 + H2} \quad (4)$$

$$h_p = SLM_r - H1 - x \quad (5)$$
<figure><img src="image_9.png" /><figcaption>**Figure 6.** (a) Comparison of cloud and shadow projections between raw and post-processed (geometrically corrected, level 1) data. (b) Geometry for the resolution of the ortho-rectification issue in two dimensions.</figcaption></figure>

Glaze et al. (1999) argued that one of the issues of the shadow technique is that the sunrays, most of the time, do not hit the cloud at the very top, but instead may intersect with the cloud a little below its top [17]. This results in an underestimation of the cloud height. In such a case, the cloud maximum height, $H_{max}$, can, however, be calculated by assuming a spherical geometry for a convective ash puff, as drawn by Sparks et al. (1997) [8]. This spherical geometry is described in Figure 7, and in such a case the height $h_s$ to add to the plume altitude of Equation (5) in order to obtain the maximum cloud height is given by:

$$h_s = 2 \sin(\frac{\alpha}{2}) \frac{d_w}{2} \cos(\frac{180 - \alpha}{2}) \qquad (6)$$

where $d_w$ is the cloud width measured with the same method as the SLM. Consequently:

$$H_{max} = \text{plume altitude (case 1 or 2)} + h_s \qquad (7)$$

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Spherical plume top geometry correction applied to estimate the maximum cloud height; a case to be used when sunrays intersect with the plume not at the very top of the plume, but a little below.</figcaption></figure>

### 3.1.4. Error Calculation

A few error sources can be neglected because of the relatively low scan angles of Low Earth Orbit (LEO) sensors such as Landsat [49]. This minimizes problems due to pixel distortion, overlap, rotation and spherical shape [48]. The first major component of the possible errors is radiation smearing. Indeed, the spectral radiance arriving at the sensor is not 100% exclusive to each pixel. The spatial distribution of radiance across the pixel is described by the point spread function [48]. This has a bell
shape, and for TM pixels spreads into one neighboring pixel [50]. The effect is especially strong when the contrasts between radiances are great, such as between dark shadow and light-colored ground. Consequently, spectral radiance from shadow pixels contaminates two to three neighboring pixels, causing the limit to blur by ±2 or 3 pixels. Moreover, the rotation performed in order to facilitate the shadow length measurements involves a restructuring of the pixels in relation to each other using cubic convolution. This operation is known to change the gray-scale value of some pixels, and can result in the insertion of "fill" pixels [51], possibly modifying the pixels selected during measurements. This effect is estimated to also affect an area of ±1 pixels. We can obtain the shadow length measurement error by multiplying the number of pixels by the pixel size previously calculated, and then the cloud height error by trigonometry using the sun azimuth angle. A summary of the resulting cloud height error estimations is given in Table 2. Cloud 2's altitude error is greater than cloud 1 because of its greater transparency, which made it more difficult to select edge and shadow pixels. Hmax is the parameter with the greatest error because its calculation used two more pixel selections (for the plume width).

<table><caption>**Table 2.** Plume height error estimations.</caption><thead><tr><th colspan="2">Cloud 1</th><th>Cloud 2</th></tr></thead><tbody><tr><td>Cloud altitude error<br/>±140 m</td><td>H<sub>max</sub> error<br/>±210 m</td><td>Cloud altitude error<br/>±185 m</td></tr></tbody></table>

We note that our solution for the effect of ortho-rectification is only for a 2D case, when the problem is three-dimensional. One solution would be the application of the Fmask algorithm developed by Zhu and Woodcock (2011) [52], and its upgrade for application to mountainous areas MFmask [40]. These algorithms are designed to automatically detect cloud and cloud shadow and match them according to their shapes. One of their inputs, though, is the cloud base height, which is needed to produce a double projection and increase the match between cloud and cloud shadow. To do this, Zhu and Woodcock (2011) used the temperature method to extract height [52]. This introduces further uncertainty inherent in applying the height-from-temperature method, which we strive to avoid here. Nevertheless, given the topography, sun angles and satellite-viewing geometry for the case considered here, we estimate only a slight underestimation of shadow length (by no more than 10%) for these low cloud heights. Error will, though, increase with cloud height, thus becoming more of a concern for higher (sub-Plinian and Plinian) clouds.

### 3.1.5. Cloud Dispersal Velocity Calculation

The ground-based video camera footage was used to calculate the cloud dispersal velocity. Frames from these data are stamped with GPS times, so that timings are accurate to less than a second. Every 10 min, two consecutive camera frames were selected (separated by less than 30 s). A point in each frame was chosen based on the same recognizable plume top shape and the number of pixels between the two points was counted. The points were chosen to be the highest and farthest from the vent possible in order to avoid the cloud dispersal velocity being affected by the jet velocity. Now velocity can be obtained from the time difference between the two frames and the distance moved. This gave a mean dispersal velocity of 18.5 m s⁻¹, with a range of 12.6 to 26.7 m s⁻¹ (Table 3). This is in agreement with the value used by De Michele et al. (2019) [38] and the wind speed profile from the National Center for Environmental Prediction (NCEP)/National Center for Atmospheric Research (NCAR) at 12:00 UTC (i.e., 18 m s⁻¹ at 10,000 m a.s.l.).

### 3.1.6. Generation of Cloud Height Time Series

Each cloud height calculated along the intersection line was linked to an emission time using the distance from the vent (d<sub>VIx</sub>) and the cloud dispersal velocity (V). To do this, the time of release of the cloud at any given point can be obtained by subtracting the time the cloud has taken to travel from
the vent to the point (=d<sub>VIX</sub>/V) from the time of image acquisition (i.e., 09:37:47 UTC). For cloud 1, the cloud dispersal velocity used was that calculated between 9:30 and 9:40. Since the proximal edge of cloud 1 was calculated to have been emitted at 9:32:58 and the image was acquired just before 09:38, only one cloud dispersal velocity was available (Table 1). For cloud 2, the cloud dispersal velocity was changed whenever the velocity for each 10-min period of cloud dispersal velocity calculation from the ECV camera changed (Table 1). The time series generated in this way from a single, static, high-spatial resolution image gives us the possibility to track plume height for several hours prior to image acquisition, with a temporal resolution of (in our case) down to 1 s over around one hour.

<table><caption>**Table 3.** Cloud dispersal velocities calculated with ECV camera footage.</caption><thead><tr><td>Time Range</td><td>Cloud Dispersal Velocity (m/s) (This Study)</td></tr></thead><tbody><tr><td>8:40–8:50</td><td>17.4</td></tr><tr><td>8:50–9:00</td><td>19.6</td></tr><tr><td>9:00–9:10</td><td>16.7</td></tr><tr><td>9:10–9:20</td><td>21.5</td></tr><tr><td>9:20–9:30</td><td>26.2</td></tr><tr><td>9:30–9:40</td><td>26.7</td></tr></tbody></table>

## 3.2. Cloud Height Validation

### 3.2.1. Cloud Height from the INGV-OE ECV Camera

The ground-based camera that we use is named after its location label (i.e., ECV) and is located 26.7 km to the South of the NSEC (see [28] for location and further camera details). Scollo et al. (2014) have developed a method to retrieve plume height from images acquired by the camera located at ECV [53], which involves placing a grid showing a succession of heights over the image. The ECV camera can capture videos of Etna's eruption columns up to a height of 9 km a.s.l, with a height calculation error of ±500 m [53]. Scollo et al. (2014) provided two examples of grids applied to two previous events [53]. These were modified in order to fit the 26 October 2013 event (Figure 8). The wind on 26 October was blowing the plume in the WSW direction while the two example grids were fitted for plumes spreading in the ESE direction. Consequently, the grid was horizontally mirrored and overlain onto the camera frames (Figure 8). The plume height was then measured as the maximum apparent height on the frame.

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Plume height retrieval from ECV camera frames (at 08:53:36 and 09:14:36) using Scollo et al.'s (2014) [53] grid (green lines). The blue lines are Etna's outline and the direction of the plume ascent (wind direction dependent). The meteorological and volcanic clouds are highlighted respectively in grey and black hatching.</figcaption></figure>
The fountaining activity started to produce a plume at around 02:45, and the sun rose around 05:00. Consequently, being a visible camera, ECV was not usable for the first two hours of the eruption. However, once enhanced, all frames after 04:30 had sufficient luminosity to allow the plume top to be viewed. As such, we were able to derive a cloud height time series with a temporal resolution of one frame per minute between 04:30 and 09:37 for validation of the OLI-derived time series.

### 3.2.2. Cloud Height from PEM

De Michele et al. (2019) created a method to extract the cloud top height from ortho-rectified Landsat 8 OLI images [38]. This method takes advantage of the physical offset between the multispectral and the panchromatic bands of the push broom sensor. This results in a short time lag (0.52 s) between acquisitions of each band. This is equivalent to the use of the stereoscopic parallax introduced by Prata and Turner (1997) [54] to take advantage of the two viewing angles provided by the Along-Track Scanning Radiometer (ATSR), but using only one single satellite pass. The Plume (top) Elevation Model (PEM) method was introduced by De Michele et al. (2016) [55] and applied to the plume from the 2014 Holuhraun (Iceland) eruption using raw (level 0) OLI data [55], and has since been developed for standard (level 1) OLI data for the same 26 October 2013 Etna lava fountain event considered here [38]. The method produces a digital elevation model for the plume top (a PEM). A new PEM for Etna's 26 October 2013 cloud was created here to correct for the ground topography. Cloud height was then retrieved at the same points used by the shadow technique. The PEM calculations, however, use a pixel window of 48 pixels and, as the window is moved across the image in steps of 24 pixels, the pixel window covers half of its neighboring window so that values are smoothed. Consequently, the spatial resolution decreases compared to the original image (with an equivalent pixel size of 240 m), and plume heights are smoothed.

## 4. Results

From here on, our cloud measurements are referred to as altitude if given in terms of their position above sea level or height if above ground level. For cloud 1, 84 height measurements were made over the 7.7 km of its downwind extent (Figure 9a). There is a gap of 22 km where no cloud is apparent. Then, cloud 2, for which we have 62 height measurements, extends 19.5 km to the image edge (Figure 9a).

### 4.1. Cloud Height Variation with Distance

Near-vent, cloud 1 ascends at a rate of 530 m of ascent per 100 m of downwind travel to reach an altitude of 7.6–8.0 km a.s.l. (4.4–4.8 km above the vent) after 1.47 km. Thereafter, ascent is less steep (90 m of ascent per 100 m of downwind travel), so that the cloud reaches a peak altitude of 10–10.4 km a.s.l. after 4.15 km. This peak altitude is equivalent to a 6.8–7.2 km high cloud. Thereafter, the cloud descends to 7.2–7.6 km a.s.l. at 7.7 km. We assume that these trends with distance reflect plume overshooting followed by adjustment (settling) to the neutral buoyancy level (cf. [41]). There is some variation in altitude after the peak at 4.15 km. For example, there are local peaks of 9–9.4 and 8.4–8.8 km at 4.8 and 6.6 km, respectively, and a low of 7.6–8 km at 5.5 km. The shape of the profile in Figure 9a is typical of a bent-over plume, where the peak point (which will be vertically over the vent in no-wind or high-MDR conditions [28,56]) has been blown 4.15 km downwind for a bent-over angle of 67°.

Cloud 2 begins 30 km away from the vent. Overall, cloud 2 ascends in altitude from 7.8–8.2 km a.s.l. at the closest point to the vent to 9.3–9.7 km a.s.l. at the image edge (59 km from the vent). It thus ascends at a rate of around 50 m per kilometer. Cloud 2 shows less frequent and less marked variations than Cloud 1, but there is a low between 33.4 and 36.2 km where the cloud altitude decreases by 700 m to a height of 7.4–7.6 km (Figure 9a).
<figure><img src="image_12.png" /><figcaption>**Figure 9.** (a) Comparison between cloud 1 and cloud 2's altitude with distance from the vent obtained using the shadow technique and PEM methods. (b) Comparison of cloud altitude time series obtained from the OLI images using the shadow technique, ECV camera and from applying Calvari et al.'s (2018) MDR versus height relation [28]. Dashed lines are extrapolated across times for which there are no data under the assumption of plume height stability under stable MDR conditions and sustained lava fountain activity.</figcaption></figure>
## 4.2. Cloud Height Time Series

Figure 9b displays the time series of plume altitude above the vent. In distance plots, the first events chronologically are furthest from the vent and the most recent events are closest to the vent. Thus, the temporal profile is a mirror of the distance profile. Empirically, MDR and plume height have been linked with increasing MDR, resulting in increasing plume height (e.g., [8,57]). This means that, assuming that the height fluctuations measured down the axis of the cloud represent variations at the source, the general increase in cloud altitude with distance down cloud 2 can be explained by a decreasing MDR with time, where the low point would be a short period of lower MDR. That is, in general, the oldest and furthest travelled part of the cloud was emitted at a higher MDR than the youngest and nearest (to the vent) part. For this case, we can trace the evolution of the plume emission and cloud dispersal back in time for 50 min prior to the acquisition of the satellite image. While the imaged portion of cloud 2 was erupted between 08:47 and 09:11, cloud 1 emission began at 09:33 and emission was continuing at the time of image acquisition at 9:37. Given an optimally placed event and low wind speeds, we can potentially derive a time series with a duration of more than 24 h. For example, the leading edge of plume with a source in the lowermost SW corner of a 180 × 185 km OLI image will, in a 3 m s⁻¹ wind blowing to the NE, take 24 h to traverse the image.

## 4.3. Comparison with PEM and ECV Camera

In Figure 9a we compare the cloud heights from the shadow with the cloud heights generated by the PEM method. There is a good agreement between these two datasets for cloud 1, except for the more proximal part where the plume begins to become transparent. There is another difference between the two height estimations at the beginning of cloud 2 for the same reason. Cloud transparency thus seems to be responsible of underestimations of between 25% and 50% in the height using the PEM method in these two regions. In addition, the PEM does not pick up the overshooting region (the highest portion of cloud 1) because the PEM method sub-samples the pixels in the image, thus smoothing the height profile. For cloud 2, there is a great deal of variation in the PEM measurements because the measurements have been made at the border of the cloud. This is a mixed pixel effect, where pixels become mixed in two senses: vertically and horizontally. In the vertical sense, the signal from a semi-transparent cloud and that from the underlying ground become mixed. In the horizontal sense, where a pixel is only partially filled by the cloud, the pixel-integrated spectral-radiance is a product of the portion of the pixel occupied by the cloud and ground occupying the remainder of the pixel (Marsh, 1981). We can also have cases of vertical and horizontal mixing (a semi-transparent cloud that partially fills a pixel). These effects result in a much greater error than for the shadow method. The PEM underestimates the height if the pixel selected is too close to the ground, or overestimates the height if the pixel selected is closer to the center of the cloud.

The cloud altitudes from the shadow technique are also in line with those of the ECV camera (Figure 9b), meaning that the plume dispersion velocities are most likely correct, and that the satellite-retrieved altitudes appear valid. Unfortunately, the weather became cloudy from around 09:20 onwards, blocking the ECV camera view of the plume top until 09:36. The satellite-retrieved cloud altitudes are, however, able to fill this gap—Figure 9b shows a generally declining plume altitude, from around 9.3–9.7 to 7.8–8.2 km a.s.l. above the vent between 08:45 and 09:10. This is a total decline of 1500 m in 25 min, for a rate of 60 m of decline per minute. There is then a rapid waxing period between 09:33 and 09:35, when the plume picked up from around 7.3–7.7 to 10–10.4 km a.s.l. This is an increase of 2600 m in 2 min, for a rate of ascent of 1300 m per minute, or around 22 m s⁻¹.

# 5. Discussion

## 5.1. Cloud Separation: Shut-Down in Activity?

Between 8 and 30 km, we see an apparent gap in the cloud in the OLI image (Figure 1). The question is: does this relate to and signify a pause in activity? Given the cloud dispersal velocity, any such
pause would have lasted for 20 min beginning at 09:13 and ending at 09:33 (Figure 9b). Cloud cover did not allow plume top height observations from the ECV camera between 09:20 and 09:36, but a lava fountain was visible from the near-vent EMOT camera at this time. The lava fountain height decreased by half between 09:12 and 09:30. Concomitantly, the plume altitude viewed by the ECV camera decreased by 500 m between 09:11 and 09:19. The difference in plume behavior is visible on Figure 8; at 9:14, the plume is much lower than at 8:53. There was then a plume altitude increase between 09:33 and 09:35, but at no point did the eruptive activity stop. This decrease in plume altitude coincides with a decrease in levels of sulfur dioxide and ash emissions (cf. Figure 4 of [58]). The short-lived (22 min long) decrease in lava fountain intensity resulted in a plume height decrease but also a decrease in the quantity of erupted material. When contrast-enhanced, there is in fact material visible in the OLI image in the air between the two clouds. However, it appears as a haze, and is not optically thick enough to produce a shadow. Thus, the gap between the two clouds, although being coincident with a decrease in activity, does not mark a shut down. The gap in the volcanic cloud is also coincident with a zone that lacks meteorological clouds (Figure 1). This coincides with the valley of the Simeto river, where elevations fall away to around 400 m a.s.l., compared with 1000 m a.s.l. in the mountains to the west. The two volcanic clouds are also light-toned, indicating that they were rich in water vapor. We thus suggest that the atmospheric conditions across the Simeto valley prevented both meteorological and water-rich (ash-poor) volcanic clouds from condensing in this area, enhancing the “disappearance” of the cloud during the decrease in activity. This apparent absence of the volcanic cloud raises concerns regarding the detection of volcanic material by satellite in this particular area, which happens to be the flight corridor for aircraft coming into the Catania Fontanarossa airport (Figure 1).

## 5.2. Settling of Cloud Height to Neutral Buoyancy Level

The height values obtained from the shadow technique and the ECV camera images are in good agreement for the proximal cloud, there being a difference of just 20–200 m (Figure 9b). Due to wind effects, the position at which the cloud reaches its maximum height is offset by 4 km to the SW (Figure 9a). Between 4 and 8 km, the cloud settles back to the level of neutral buoyancy, which is around 1 km below the overshoot height. For the case considered here, settling is systematic and follows the trend:

$$ \text{Cloud altitude (m)} = -0.415 * \text{Distance from the vent (m)} + 11.093. (r^2 = 0.68) \quad (8) $$

It is important to note that this relation lumps all of the complexities in the plume settling dynamics into a single, simple relation, where our objective is simply to derive a trend for plume height with distance for this case where a simple, systematic trend can be seen.

Calvari et al. (2018) proposed a formula linking lava fountain, *H* (lava fountain), and associated plume, *H* (plume), heights (in m) [28]:

$$ H(\text{plume}) = 5.26 H(\text{lava fountain}) + 6.830 \quad (9) $$

This relation was empirically derived by considering data from 20 lava fountaining events at the NSEC between 2011 and 2013, but did not include the event of the case study given here. The values calculated with Equation (8) are systematically higher than the heights obtained by the shadow method (Figure 9b). However, our measurements are mostly for the height of neutral buoyancy and not maximum plume height. Thus, the 800–2200 m difference in Figure 9b between the result of Equation (9) and the measurement of plume height from the shadow may be due to plume settling. This therefore likely gives the difference between maximum plume height, *H* (plume), and neutral buoyancy height, *H* (NB), so that:

$$ H(\text{NB}) = H(\text{plume}) - \text{average offset} = 5.26 H(\text{lava fountain}) + 5.130. (r^2 = 0.50) \quad (10) $$
Here, the average offset is the median value for the differences between Equation (9)'s output and the plume height at all down-wind measurement points. This relation is supported by the fact that the maximum plume height measured from the image is, in fact, in very good agreement with that expected from the relation of Calvari et al. (2018) [28] (Figure 9b). Thereafter, the results of Equation (9) and the measurement of cloud height diverge downwind due to plume settling.

### 5.3. Mass Flux Fluctuation Effect

Assuming constant rates of air entrainment (cf. [59]), an empirical relation that enables the estimation of the MDR required to drive a plume to height *H* (*plume*) is that of Woodhouse et al. (2013) [60]:

$$H(\text{plume}) = 0.318 * \text{MDR}^{0.253} \qquad (11)$$

where *H* (*plume*) is measured in kilometers and MDR is the mass flux measured in kilograms per second. However, this relation is designed for sustained plumes associated with Plinian eruptions [60], and not for lava fountains, which have different convective properties. During lava fountaining events, magmas can ascend the conduit to erupt as a spray of magma clots and gas at the surface [8,61]. A variable—but usually minor when compared with Plinian plumes—amount of fine pyroclasts is produced to drive a convective plume above the fountain [8]. Most of the heat produced during lava fountains is held by the coarse ejecta, which fall from the plume rapidly [8]. As a result, only a small percentage of the thermal energy flux is conveyed to the plume, so that volcanic plumes above fountaining vents are usually weaker and less high than those fed at comparable MDR during Vulcanian and sub-Plinian eruptions [8]. Moreover, Equation (11) requires maximum plume height and not the neutral buoyancy height calculated in most of the points in this study. We thus need to adapt Equation (11) for a jet-fed thermal, using the height corrected for the settling effect following Equation (10). Calvari et al. (2018) provide a time series of volumetric discharge rates (VDR) for this event [28]. We used the relation of Woodhouse et al. (2013) [60] to estimate the MDR from each of our *H* (*plume*) and corrected the height using Equation (11), then converted them to VDR using an assumed tephra density of 2650 kg m⁻³ dense rock equivalent (DRE) [28] and compared them (Figure 10).

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Comparison between cloud volumetric discharge rate (VDR) estimated from *H* (*plume*) using Equation (10) from Woodhouse et al. (2013) [60] and measurements of Calvari et al. (2018) [28].</figcaption></figure>

Following Andronico et al. (2018) [37], only 4% of the erupted magma volume DRE accounts for the distal tephra fallout (i.e., the cloud), and 23% for the proximal tephra fallout (the lava fountain). This means that the mass discharge rate for particles going into the cloud should be around 17% of that for those involved in feeding the lava fountain. Considering that the VDR calculated from Calvari et al. (2018) fluctuates between 100 and 150 m³ s⁻¹ [28], the DRE flux (for the particulate portion of
the mixture) should be around 17–26 m³ s⁻¹. Instead, in Figure 10, the DRE flux fluctuates between 250 and 450 m³ s⁻¹. However, the difference is systematic, allowing us to multiply the relation of Woodhouse et al. (2013) [60] by a factor of 0.065 ± 0.005 to take into account the differing dynamics and mass partitioning between a Plinian and a lava fountain plume. The following corrected equation was only fitted to our single event, and applies to the volume flux of solid particles (as DRE) and not the whole cloud (mixture of particles, gas and entrained air):

$$H(\text{plume}) = 0.318 * (VDR * 0.065 \pm 0.005)^{0.253} \quad (12)$$

## 5.4. Wind Effects on the Plume above the Vent

The effects of the wind on the plume are visible in the cloud 1 profile (Figure 9), because the crosswind speed is higher than the plume rise velocity [25]. At the acquisition time, the plume dispersal velocity was 12.5 m s⁻¹. However, there seem to be two dynamics within the ascending column (Figure 11). It first ascends with a slope of 530 m per 100 m of downwind travel to reach an altitude of 7.6–8 km a.s.l. (4.4–4.8 km above the vent) for a downwind bending 1.47 km. Then, the ascent is less steep (90 m per 100 m of downwind travel) until the cloud reaches its peak altitude of 10–10.4 km a.s.l. (7840 m above the vent) after 4.15 km. The first part of the ascending column has an angle of 72° and the second part an angle of 51° (Figure 11).

<figure><img src="image_14.png" /><figcaption>Figure 11. Wind effect on the eruptive column of Etna 26 October 2013 lava fountain.</figcaption></figure>

The plume is thus more bent in its upper part than in its lower part (Figure 11). This fits with a non-uniform wind field with height [56] as well as a lower mass of solids feeding the upper part. In fact, the wind profile versus altitude from the NCEP/NCAR data for Etna at 12:00 UTC [38] shows that the wind speed increased with altitude between 3 km (the altitude of the vent) up to 10 km a.s.l. Thus the difference in the bending angle is the result of the transition between the fountain part of the plume, which is jet-dominated with a relatively high MDR and relatively low wind speed (Figure 11), and the upper tephra-rich part where the wind speed is higher, where the MDR is lower, and where more air is incorporated inside the plume [8,56]. The upper part has transitioned to a buoyant phase, where the air entrainment rates can be significantly higher for buoyant thermals fed by a jet in such “weak” (Strombolian and Hawaiian) plume cases [22].

## 5.5. Plume/Cloud Ascent and Dispersion Dynamics: Summary

The 26 October 2013 fountain event at Etna produced a plume that, above the vent, was a bent-over rooted thermal. A decrease in mass flux between the lower fountain-fed portion of the plume and the upper portion, in tandem with the increasing wind speed with altitude, resulted in an increase in the
degree of bending 2.5 km above the vent. This represents the transition between the jet of large clasts that fall back to the ground (thereby removing mass and thermal energy from the system) and the plume of tephra that continues to ascend by convection (Figure 11). In such a scenario, the upper part is bent over more easily than the lower part. Ascent rates were 22 m s⁻¹ up to a maximum height of 10–10.4 km above (but offset by 4.2 km downwind from) the vent. Thereafter, overshooting tephra settled to a neutral buoyancy level at around 7.8–8.2 km a.s.l. after 7.7 km (Figure 12). A lack of condensation in the water-rich plume meant that there was an apparent gap in the cloud between 8 and 30 km. However, thereafter an apparent increase in cloud height implied higher mass fluxes in the 50-min period prior to image acquisition (Figure 12).

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Summary of the dispersion dynamics of the 26 October 2013 Etna fountain-fed plume and cloud system. The red dot indicates the visible end of cloud 2 at the image edge.</figcaption></figure>

**6. Conclusions**

The cloud-height-from-shadow technique set up and validated in this study allowed the retrieval of a cloud altitude time-series from a single LANDSAT 8 OLI image, allowing us to document the ascent and dispersion history of a plume-cloud system emitted during a fountaining event at Etna volcano. The high-spatial resolution of the LANDSAT 8 product allowed us to detail and quantify cloud and plume dynamics over a distance of 60 km and over a time period of 50 min, with the precision of a few seconds and vertical error on plume altitude of ±200 m. Potentially, our method allows downwind plume height profiles and MDR time series to be built over distances of up to 260 km and periods of 24 h, depending on source location in the image, wind speed and direction. The results were found to be in good agreement with measurements of plume height from ground-based cameras, as well as the PEM method of De Michele et al. (2019) [38].

The data set derived here for cloud/plume altitude versus distance and time allowed a number of empirical relations to be derived and refined for a fountain-fed plume. These include refinement for relations between settling and distance, maximum plume height and level of neutral buoyancy, and MDR and plume height. These relations are, however, strictly calibrated for the 26 October 2013 event, but do suggest peculiarities in the dynamics and mass partitioning between different segments of a fountain-fed plume when compared with Vulcanian and Plinian plumes. However, whether the empirical, best-fitting relations given here can be viewed as a general model for fountain-fed plumes requires the collection of data for further events so as to build a model based on a robust statistical
approach. We thus advocate the targeting of high-spatial resolution satellite sensors to fountaining events, especially those where ancillary ground-based thermal camera videos are available (e.g., Etna, Hawaii, and Piton de la Fournaise) to follow up on, and check, the relations implied here.

The time series generated here uses a single, static, high-spatial resolution image, allowing time series with a resolution of seconds to be built from a single image. The temporal resolution depends on the distance between pixels selected, whereby an increase in the number of points where measurements are made will increase the processing time. Potentially, though, the measurement can be made pixel-by-pixel, allowing a measurement to be made every few tens of meters which, in high wind speed conditions, will convert to a measurement every few seconds. Such time series have, to date, typically been the domain of sensors mounted on geostationary satellites with spatial resolutions of 1–4 km, and temporal resolutions of 15 min (e.g., [16,21,42]). The method thus represents a great increase in both the spatial and temporal resolutions of the record, providing data essential for constraining and running models for plume ascent and dispersal cf. [62]. There are, however, a few limits to the application of this method. First, it can only apply to daytime events, as we required the presence of a shadow. Second, the conversion to a time series requires reliable data on wind speed or cloud dispersal velocity. Because these vary with time, height and distance, a single value might not give reliable results. Finally, such explosive events (being short, just a few hours in duration) are difficult to capture given the typical 16-day return period of LEO satellites carrying high spatial resolution sensors needed for the method applied here. This argues for more use of pointing capabilities, such as that represented by the ASTER Urgent Response Protocol [63], or constellations, wherein currently the Landsat-ASTER-Sentinel satellites represent just such a network.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2072-4292/12/23/3951/s1, Table S1: Reference list for volcanological applications of satellite data to volcanic plumes. Table S2: Sun-Earth geometry calculations.

**Author Contributions:** Conceptualization, A.J.L.H. and L.G.; methodology, S.P.-B., A.J.L.H., S.C. and M.D.M.; validation, A.J.L.H. and S.C.; formal analysis, S.P.-B., S.C. and M.D.M.; investigation, S.P.-B., A.J.L.H. and S.C.; resources, S.P.-B., A.J.L.H., S.C. and, M.D.M.; data curation, S.P.-B.; writing—original draft preparation, S.P.-B.; writing—review and editing, A.J.L.H., S.C., M.D.M. and L.G.; visualization, S.P.-B.; supervision, A.J.L.H.; project administration, L.G.; funding acquisition, L.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by CNES-TOSCA (Terre Solide), grant number 10 3703 "Integration of sample return data and remote sensing for advanced understanding of volcanic ash formation and dispersion" (PI: Lucia Gurioli).

**Acknowledgments:** The authors would like to thank the INGV-OE personnel for the cameras network maintenance. This manuscript benefited greatly from the valuable input on satellite viewing geometry and orthorectification of Daniel Raucoules (BRGM Orléans). We are grateful to Mike Ramsey and Diego Coppola who assisted us in satellite data collection. We owe much to Nathan Roberts from USGS support staff and the "calibration/validation team" for details on Landsat's acquisition time information. We wish to thank Stefano Mannini for his input on satellite data processing.

**Conflicts of Interest:** The authors declare no conflict of interest.

## References

1. Casadevall, T.J. *Volcanic Ash and Aviation Safety*; Proceedings of the First International Symposium on Volcanic Ash and Aviation Safety; U.S. Geological Survey Bulletin: Washington, DC, USA, 1994.

2. Kienle, J.; Shaw, G.E. Plume dynamics, thermal energy and long-distance transport of vulcanian eruption clouds from Augustine Volcano, Alaska. *J. Volcanol. Geotherm. Res.* **1979**, *6*, 139–164. [CrossRef]

3. Hanstrum, B.N.; Watson, A.S. A case study of two eruptions of Mount Galunggung and an investigation of volcanic eruption cloud characteristics using remote sensing techniques. *Australian Meteorol. Mag.* **1983**, *31*, 171–177.

4. Aloisi, M.; D'Agostino, M.; Dean, K.G.; Mostaccio, A.; Neri, G. Satellite analysis and PUFF simulation of the eruptive cloud generated by the Mount Etna paroxysm of 22 July 1998. *J. Geophys. Res. Solid Earth* **2002**, *107*, ECV 9-1–ECV 9-12. [CrossRef]
5. Guffanti, M.; Mayberry, G.C.; Casadevall, T.J.; Wunderman, R. Volcanic hazards to airports. *Nat. Hazards* **2009**, *51*, 287-302. [CrossRef]

6. Alemanno, A. *Governing Disasters—The Challenges of Emergency Risk Regulation*; Edward Elgar Publishing Limited: Cheltenham, UK, 2011; ISBN 9780857935731.

7. Oppenheimer, C. Volcanological applications of meteorological satellites. *Int. J. Remote Sens.* **1998**, *19*, 2829-2864. [CrossRef]

8. Sparks, R.S.J.; Bursik, M.; Carey, S.; Gilbert, J.S.; Glaze, L.S.; Sigurdsson, H.; Woods, A.W. *Volcanic Plumes*; Wiley: Chichester, UK, 1997; ISBN 9780471939016.

9. Bonadonna, C.; Folch, A.; Loughlin, S.; Puempel, H. Future developments in modelling and monitoring of volcanic ash clouds: Outcomes from the first IAVCEI-WMO workshop on Ash Dispersal Forecast and Civil Aviation. *Bull. Volcanol.* **2012**, *74*, 1-10. [CrossRef]

10. Newhall, C.G.; Self, S. The volcanic explosivity index (VEI) an estimate of explosive magnitude for historical volcanism. *J. Geophys. Res.* **1982**, *87*, 1231. [CrossRef]

11. Francis, P.W. Infra-red techniques for volcano monitoring and prediction-A review. *J. Geol. Soc. Lond.* **1979**, *136*, 355-359. [CrossRef]

12. Francis, P.W.; Oppenheimer, C. Applications of satellite remote sensing techniques to volcanology. In *Understanding the Terrestrial Environment—The Role of Observation from Space*; Paylor & Francis: London, UK, 1992; pp. 37-52.

13. Francis, P.W.; Oppenheimer, C. *Volcanoes*, 2nd ed.; Oxford University Press: Oxford, UK, 2003; ISBN 9780199254699.

14. Pyle, D.M. Sizes of Volcanic Eruptions. In *The Encyclopedia of Volcanoes*; Elsevier: London, UK, 2015; pp. 257-264.

15. Mouginis-Mark, P.J.; Pieri, D.C.; Francis, P.W.; Wilson, L.; Self, S.; Rose, W.I.; Wood, C.A. Remote sensing of volcanos and volcanic terrains. *EOS Trans. Am. Geophys. Union* **1989**, *70*, 1567-1575. [CrossRef]

16. Glaze, L.S.; Francis, P.W.; Self, S.; Rothery, D.A. The 16 September 1986 eruption of Lascar volcano, north Chile: Satellite investigations. *Bull. Volcanol.* **1989**, *51*, 149-160. [CrossRef]

17. Glaze, L.S.; Wilson, L.; Mouginis-Mark, P.J. Volcanic eruption plume topography and heights as determined from photoclinometric analysis of satellite data. *J. Geophys. Res. Solid Earth* **1999**, *104*, 2989-3001. [CrossRef]

18. Woods, A.W.; Self, S. Thermal disequilibrium at the top of volcanic clouds and its effect on estimates of the column height. *Nature* **1992**, *355*, 628-630. [CrossRef]

19. Holasek, R.E.; Self, S.; Woods, A.W. Satellite observations and interpretation of the 1991 Mount Pinatubo eruption plumes. *J. Geophys. Res. B Solid Earth* **1996**, *101*, 27635-27655. [CrossRef]

20. Zakšek, K.; Hort, M.; Zaletelj, J.; Langmann, B. Monitoring volcanic ash cloud top height through simultaneous retrieval of optical data from polar orbiting and geostationary satellites. *Atmos. Chem. Phys.* **2013**, *13*, 2589-2606. [CrossRef]

21. Marchese, F.; Falconieri, A.; Pergola, N.; Tramutoli, V. A retrospective analysis of the Shinmoedake (Japan) eruption of 26-27 January 2011 by means of Japanese geostationary satellite data. *J. Volcanol. Geotherm. Res.* **2014**, *269*, 1-13. [CrossRef]

22. Patrick, M.R. Dynamics of Strombolian ash plumes from thermal video: Motion, morphology, and air entrainment. *J. Geophys. Res.* **2007**, *112*, B06202. [CrossRef]

23. Carey, S.; Bursik, M. Volcanic Plumes. In *The Encyclopedia of Volcanoes*; Elsevier: London, UK, 2015; pp. 571-585.

24. Bonaccorso, A.; Calvari, S. A new approach to investigate an eruptive paroxysmal sequence using camera and strainmeter networks: Lessons from the 3-5 December 2015 activity at Etna volcano. *Earth Planet. Sci. Lett.* **2017**, *475*, 231-241. [CrossRef]

25. Carey, S.; Sparks, R.S.J. Quantitative models of the fallout and dispersal of tephra from volcanic eruption columns. *Bull. Volcanol.* **1986**, *48*, 109-125. [CrossRef]

26. Slawson, P.R.; Csanady, G.T. On the mean path of buoyant, bent-over chimney plumes. *J. Fluid Mech.* **1967**, *28*, 311. [CrossRef]

27. Andò, B.; Pecora, E. An advanced video-based system for monitoring active volcanoes. *Comput. Geosci.* **2006**, *32*, 85-91. [CrossRef]
28. Calvari, S.; Cannavò, F.; Bonaccorso, A.; Spampinato, L.; Pellegrino, A.G. Paroxysmal Explosions, Lava Fountains and Ash Plumes at Etna Volcano: Eruptive Processes and Hazard Implications. *Front. Earth Sci.* **2018**, *6*. [CrossRef]

29. Behncke, B.; Branca, S.; Corsaro, R.A.; De Beni, E.; Miraglia, L.; Proietti, C. The 2011–2012 summit activity of Mount Etna: Birth, growth and products of the new SE crater. *J. Volcanol. Geotherm. Res.* **2014**, *270*, 10–21.
[CrossRef]

30. De Beni, E.; Behncke, B.; Branca, S.; Nicolosi, I.; Carluccio, R.; D'Ajello Caracciolo, F.; Chiappini, M. The continuing story of Etna's New Southeast Crater (2012–2014): Evolution and volume calculations based on field surveys and aerophotogrammetry. *J. Volcanol. Geotherm. Res.* **2015**, *303*, 175–186. [CrossRef]

31. Coltelli, M.; Del Carlo, P.; Vezzoli, L. Discovery of a Plinian basaltic eruption of Roman age at Etna volcano, Italy. *Geology* **1998**, *26*, 1095. [CrossRef]

32. Istituto Nazionale di Statistica-Roma, I. Demografia in Cifre. Available online: http://demo.istat.it/ (accessed on 1 October 2020).

33. Vulpiani, G.; Ripepe, M.; Valade, S. Mass discharge rate retrieval combining weather radar and thermal camera observations. *J. Geophys. Res. Solid Earth* **2016**, *121*, 1–17. [CrossRef]

34. Del Carlo, P.; Vezzoli, L.; Coltelli, M. Last 100 ka tephrostratigraphic record of Mount Etna. *GMS* **2004**, *143*, 77–89.

35. Corradini, S.; Guerrieri, L.; Lombardo, V.; Merucci, L.; Musacchio, M.; Prestifilippo, M.; Scollo, S.; Silvestri, M.; Spata, G.; Stelitano, D. Proximal Monitoring of the 2011–2015 Etna Lava Fountains Using MSG-SEVIRI Data. *Geosciences* **2018**, *8*, 140. [CrossRef]

36. Greco, F.; Currenti, G.; Palano, M.; Pepe, A.; Pepe, S. Evidence of a shallow persistent magmatic reservoir from joint inversion of gravity and ground deformation data: The 25–26 October 2013 Etna lava fountaining event. *Geophys. Res. Lett.* **2016**, *43*, 3246–3253. [CrossRef]

37. Andronico, D.; Behncke, B.; De Beni, E.; Cristaldi, A.; Scollo, S.; Lopez, M.; Lo Castro, M.D. Magma Budget From Lava and Tephra Volumes Erupted during the 25-26 October 2013 Lava Fountain at Mt Etna. *Front. Earth Sci.* **2018**, *6*. [CrossRef]

38. De Michele, M.; Raucoules, D.; Corradini, S.; Merucci, L.; Salerno, G.; Sellitto, P.; Carboni, E. Volcanic cloud top height estimation using the plume elevation model procedure applied to orthorectified Landsat 8 data. test case: 26 October 2013 Mt. Etna eruption. *Remote Sens. 2019*, *11*, 785. [CrossRef]

39. Young, N.E.; Anderson, R.S.; Chignell, S.M.; Vorster, A.G.; Lawrence, R.; Evangelista, P.H. A survival guide to Landsat preprocessing. *Ecology* **2017**, *98*, 920–932. [CrossRef]

40. Qiu, S.; He, B.; Zhu, Z.; Liao, Z.; Quan, X. Improving Fmask cloud and cloud shadow detection in mountainous area for Landsats 4–8 images. *Remote Sens. Environ.* **2017**, *199*, 107–119. [CrossRef]

41. Dean, K.G.; Bowling, S.A.; Shaw, G.E.; Tanaka, H. Satellite analyses of movement and characteristics of the Redoubt Volcano plume, January 8, 1990. *J. Volcanol. Geotherm. Res.* **1994**, *62*, 339–352. [CrossRef]

42. Holasek, R.E.; Self, S. GOES weather satellite observations and measurements of the May 18, 1980, Mount St. Helens eruption. *J. Geophys. Res. Solid Earth* **1995**, *100*, 8469–8487. [CrossRef]

43. Kinoshita, K. Observation of flow and dispersion of volcanic clouds from Mt. Sakurajima. *Atmos. Environ.* **1996**, *30*, 2831–2837. [CrossRef]

44. Denniss, A.M.; Harris, A.J.L.; Rothery, D.A.; Francis, P.W.; Carlton, R.W.T. Satellite observations of the April 1993 eruption of Lascar volcano. *Int. J. Remote Sens.* **1998**, *19*, 801–821. [CrossRef]

45. Prata, A.J.; Grant, I.F. Retrieval of microphysical and morphological properties of volcanic ash plumes from satellite data: Application to Mt Ruapehu, New Zealand. *Q. J. R. Meteorol. Soc.* **2001**, *127*, 2153–2179.
[CrossRef]

46. Tupper, A.; Carn, S.A.; Davey, J.; Kamada, Y.; Potts, R.; Prata, F.; Tokuno, M. An evaluation of volcanic cloud detection techniques during recent significant eruptions in the western 'Ring of Fire'. *Remote Sens. Environ.* **2004**, *91*, 27–46. [CrossRef]

47. Iqbal, M. Sun-Earth Astronomical Relationships. In *An Introduction to Solar Radiation*; Elsevier: London, UK, 1983; pp. 1–28. ISBN 9780123737502.

48. Harris, A.J.L. *Thermal Remote Sensing of Active Volcanoes*; Cambridge University Press: Cambridge, UK, 2013; ISBN 9781139029346.

49. Mouginis-Mark, P.J.; Domergue-Schmidt, N. Acquisition of satellite data for volcano studies. *Geophys. Monogr. Ser.* **2000**, *116*, 9–24. [CrossRef]
50. Markham, B.L. Characterization of the Landsat sensors' spatial responses. *IEEE Trans. Geosci. Remote Sens.* **1985**, *GE-23*, 864-875. [CrossRef]

51. Oppenheimer, C.; Francis, P.W.; Rothery, D.A.; Carlton, R.W.T.; Glaze, L.S. Infrared image analysis of volcanic thermal features: Lascar Volcano, Chile, 1984-1992. *J. Geophys. Res.* **1993**, *98*, 4269-4286. [CrossRef]

52. Zhu, Z.; Woodcock, C.E. Object-based cloud and cloud shadow detection in Landsat imagery. *Remote Sens. Environ.* **2012**, *118*, 83-94. [CrossRef]

53. Scollo, S.; Prestifilippo, M.; Pecora, E.; Corradini, S.; Merucci, L.; Spata, G.; Coltelli, M. Eruption column height estimation of the 2011-2013 Etna lava fountains. *Ann. Geophys.* **2014**, *57*, [CrossRef]

54. Prata, A.J.; Turner, P.J. Cloud-top height determination using ATSR data. *Remote Sens. Environ.* **1997**, *59*, 1-13. [CrossRef]

55. De Michele, M.; Raucoules, D.; Arason, Þ. Volcanic Plume Elevation Model and its velocity derived from Landsat 8. *Remote Sens. Environ.* **2016**, *176*, 219-224. [CrossRef]

56. Bursik, M. Effect of wind on the rise height of volcanic plumes. *Geophys. Res. Lett.* **2001**, *28*, 3621-3624. [CrossRef]

57. Mastin, L.G.; Guffanti, M.; Servranckx, R.; Webley, P.W.; Barsotti, S.; Dean, K.G.; Durant, A.; Ewert, J.W.; Neri, A.; Rose, W.I.; et al. A multidisciplinary effort to assign realistic source parameters to models of volcanic ash-cloud transport and dispersion during eruptions. *J. Volcanol. Geotherm. Res.* **2009**, *186*, 10-21. [CrossRef]

58. Sellitto, P.; Di Sarra, A.; Corradini, S.; Boichu, M.; Herbin, H.; Dubuisson, P.; Sèze, G.; Meloni, D.; Monteleone, F.; Merucci, L.; et al. Synergistic use of Lagrangian dispersion and radiative transfer modelling with satellite and surface remote sensing measurements for the investigation of volcanic plumes: The Mount Etna eruption of 25-27 October 2013. *Atmos. Chem. Phys.* **2016**, *16*, 6841-6861. [CrossRef]

59. Suzuki, Y.J.; Koyaguchi, T. Effects of wind on entrainment efficiency in volcanic plumes. *J. Geophys. Res. Solid Earth* **2015**, *120*, 6122-6140. [CrossRef]

60. Woodhouse, M.J.; Hogg, A.J.; Phillips, J.C.; Sparks, R.S.J. Interaction between volcanic plumes and wind during the 2010 Eyjafjallajökull eruption, Iceland. *J. Geophys. Res. Solid Earth* **2013**, *118*, 92-109. [CrossRef]

61. Allard, P.; Burton, M.; Muré, F. Spectroscopic evidence for a lava fountain driven by previously accumulated magmatic gas. *Nature* **2005**, *433*, 407-410. [CrossRef]

62. Costa, A.; Suzuki, Y.J.; Cerminara, M.; Devenish, B.J.; Ongaro, T.E.; Herzog, M.; Van Eaton, A.R.; Denby, L.C.; Bursik, M.; de' Michieli Vitturi, M.; et al. Results of the eruptive column model inter-comparison study. *J. Volcanol. Geotherm. Res.* **2016**, *326*, 2-25. [CrossRef]

63. Ramsey, M.S. Synergistic use of satellite thermal detection and science: A decadal perspective using ASTER. *Geol. Soc. Lond. Spec. Publ.* **2016**, *426*, 115-136. [CrossRef]

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_16.png" /></figure>

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).