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

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.


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][2][3][4][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

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

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.

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 purposed 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.

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).

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][42][43][44][45][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.

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.
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.

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: where α is the sun elevation and σ 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

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: where α is the sun elevation and σ 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 Remote Sens. 2020, 12, 3951 7 of 23 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.
Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 23 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.

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: ( 1) = tan 1 tan ( sin ( 1 )) + ( 2) = tan 1 tan ( sin ( 1 )) + 1 tan + 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.

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: 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.
Remote Sens. 2020, 12, 3951 9 of 23 Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 23 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 (SLMr), i.e., corrected for geometric distortion, if we are to correctly calculate the actual cloud height, hp ( 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 SLMr (Equation (4)), which then allows us to estimate hp using Equation (5): 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 (SLM 1 ) data, but we need the shadow length measurement in reality (SLM r ), i.e., corrected for geometric distortion, if we are to correctly calculate the actual cloud height, h p (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 r (Equation (4)), which then allows us to estimate h p using Equation (5): 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, Hmax, 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 hs to add to the plume altitude of Equation (5) in order to obtain the maximum cloud height is given by: where dw is the cloud width measured with the same method as the SLM. Consequently: 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.

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  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: where d w is the cloud width measured with the same method as the SLM. Consequently: Remote Sens. 2020, 12, x FOR PEER REVIEW 10 of 23 (a) (b) 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, Hmax, 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 hs to add to the plume altitude of Equation (5) in order to obtain the maximum cloud height is given by: where dw is the cloud width measured with the same method as the SLM. Consequently: 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.

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

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. H max is the parameter with the greatest error because its calculation used two more pixel selections (for the plume width). 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.

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 −1 , with a range of 12.6 to 26.7 m s −1 (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 −1 at 10,000 m a.s.l.).

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 VIx ) 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 VIx /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.

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.

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 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.

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.

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).

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 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 -1 wind blowing to the NE, take 24 h to traverse the image.

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 −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).

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: Cloud altitude (m) = −0.415 * Distance f rom the vent (m) + 11.093. r 2 = 0.68 (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]: 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: 26 H(lava f ountain) + 5.130. r 2 = 0.50 (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.

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]: 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 −3 dense rock equivalent (DRE) [28] and compared them ( Figure 10).  (Figure 9b). Thereafter, the results of Equation (9) and the measurement of cloud height diverge downwind due to plume settling.

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]: 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 −3 dense rock equivalent (DRE) [28] and compared them ( Figure  10). 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 3 s −1 [28], the DRE flux (for the particulate portion of the mixture) should be around 17-26 m 3 s −1 . Instead, in Figure 10 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 3 s −1 [28], the DRE flux (for the particulate portion of the mixture) should be around 17-26 m 3 s −1 . Instead, in Figure 10, the DRE flux fluctuates between 250 and 450 m 3 s −1 . 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(plume) = 0.318 * (VDR * 0.065 ± 0.005) 0.253 (12)

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 −1 . 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.

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 −1 . 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.  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].

Plume/Cloud Ascent and Dispersion Dynamics: Summary
The 26 October 2013 fountain event at Etna produced a plume that, above the vent, was a bentover 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 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].

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 −1 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).
Remote Sens. 2020, 12, x FOR PEER REVIEW  19 of 23 scenario, the upper part is bent over more easily than the lower part. Ascent rates were 22 m s −1 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).

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.

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.