*3.2. Remote Sensing Data (EO-1 (ALI) and SPOT 6*/*7) Used for the Interpretation of the Eruptive Temporal Sequence*

SPOT6/7 and EO-1 (ALI) images were used to support the chronology of the eruptive phase during the period of 2014–2016 (Table 2). The technical characteristics of the SPOT 6/7 satellites were mentioned above. ALI is one of the instruments onboard the EO-1 satellite. This is a multispectral sensor of 10 bands: R, G, B, VNIR, SWIR (30 m of resolution) and Panchromatic (10 m resolution), the satellite provides a standard swath around 37 km and length of 42 km (and it is possible to get increased length of 185 km); and a revisit on ground each 26 days. The EO-1 was part of NASA's project "Volcano Sensor Web" [33], focused on monitoring volcanic activity in selected areas including the Kilauea volcano, Iceland's Eyjafjallajökull and Anak Krakatau Volcano; EO-1 was used principally for tracking lava flows, discharge-rate estimation relative and thermal remote sensing for hotspot detection [33–36].


**Table 2.** Remote Sensing and SPOT data used for the eruptive temporal sequence.

All SPOT images were acquired in L1-A level, EO-1 (ALI) image in L1-B format. Radiometric corrections were applied using the "Apply Gain and Offset" algorithm and FLAASH atmospheric correction module of the ENVI 5.3 software. These corrections consist of removing atmospheric distortions due to the presence of water vapor and aerosols based on estimations of the correct wavelength positions of the bands. To achieve an accurate coregistration of each DSM, a 5-m DSM based on light detection and ranging data (LIDAR) was used as a spatial reference. The LIDAR data were acquired in 2012 (2012-DSM) and obtained from the Mexican National Institute of Statistics and Geography [37]. Moreover, to improve the visual interpretation of the images, spectral enhancement was applied to maximize the contrast of the data by applying a non-linear contrast stretch-enhanced RGB composite image [38]. Specifically, from the 2015 SPOT image, a 2.5 m resolution image was obtained with a fusion process between the 10 m resolution multispectral image and the panchromatic data (Figure 2b), which, in combination with the E0-1 (ALI) data, was used to better define the spatial distribution of the lava flows prior to and after the climatic activity of July 2015 (Figure 2).

#### **4. Methods**

The work flow applied for the volume estimation of the lava and pyroclastic flows is shown in Figure 4, which is described in the sub-sections below.

**Figure 4.** Work flow applied for volume estimations of the lava and pyroclastic flows linked to the eruptive phase (2014–2016). DSM-LIDAR, digital surface model light detection and ranging data; WLF, lava flow on the W flank; SWLF, lava flow on the SW slope; SLF, south lava flow.

#### *4.1. Digital Surface Models (DSM) Processing*

Toperform thedigital stereoscopyprocessing,weused the IMAGINEPhotogrammetry 2015 software. This is a specialized module of ERDAS IMAGINE that focuses on high-accuracy full photogrammetry processing, such as digital terrain model generation, full analytical triangulation and orthophoto-mosaicking production. We used an intuitive interface with a linear workflow that was based on the model sensors that were previously defined for each aircraft and satellite sensor. In accordance with the SPOT image, the rational polynomial coefficients (RPC) file contains all of the orbital ephemeris coefficients required to resolve the software model and thus acquire the point cloud needed to generate the final DSM. To improve the absolute horizontal (x, y) position, an ortho-rectified photo mosaic with a pixel resolution of 2 m was used as a reference. This mosaic was acquired from the National Institute of Geography and Statistics (INEGI) and was, therefore, the highest quality georeference in Mexico. Approximately 25 and 30 ground control points (GCPs) were settled for each generated DSM. Twenty tie points were also selected to increase the spatial correlation between the stereo and tri-stereo images. Finally, three DSMs were extracted using the WGS 84 datum and the Universal Transversal of Mercator (UTM) projection with a horizontal resolution between 7 m and 21 m (Table 1). From this point on, the three DSMs will be abbreviated as follows: 2014-DSM, 2015-DSM, and 2017-DSM.

To evaluate the DSM accuracies, the root mean squared error (RMSE) was calculated for each DSM. The highest National Geodesic Network accuracy was used as a spatial reference; this network is a collection constituted by more than 100,000 geodetic control points distributed over all territory since 1978 (Table 3). The 2015- and 2017-DSMs have RMSEs below 1 m; in contrast, the 2014-DSM extracted from the tri-stereo 2014 data has a RMSE of 1.29 m, which represents the lowest accuracy among the three DSMs extracted with the stereoscopic techniques. We would have expected a higher vertical accuracy of the DSMs extracted from the tri-stereo images due to the tri-stereoscopic capacity, which would reduce the number of voids caused by the high topography in volcanic areas. However, this result can be explained by the small angle formed between the triplet image acquisitions (+19.08◦/+8.22◦/+14.05◦) with respect to an optimal tri-stereo configuration of −15◦–20◦/±0◦/+15◦+20◦.


**Table 3.** Root mean squared error (RMSE) estimation based on the National Geodesic Network of Mexico.

### *4.2. Volume Estimation of the Associated Lava Flows and Pyroclastic Flows*

As mention before, a 5 m DSM from the LIDAR data acquired in 2012 (2012-DSM) was used as a pre-eruptive DSM topographic reference. For the estimation of the lava flow volumes, the 2012-DSM and 2017-DSM were used as pre- and post-topography references to obtain the volumes of the main lava flows that were emplaced on the cone during the 2014–1016 period.

The estimation was performed with raster operations using @ArcGis 10.2 (Table 4, Figure 5). The 2014-DSM was used only to evaluate the partial advance of the SW lava flow that, at the time of the 2014 SPOT image, was still advancing downwards (Figure 6a,b). Specifically, for the two lava flows that were emplaced in May and September 2014, the 2012-DSM was appropriated as the pre-event surface because, between 2012 and 2014, no major events affected the W or SW sectors of the cone besides the emplacement of small pyroclastic flows and lavas (Figure 5, black arrow). In contrast, for the lava flow that was emplaced in September 2016 on the S flank, the 2015-DSM could not be used as a pre-event surface because the cone was covered by clouds. Thus, the 2012-DSM was used, but this was still inappropriate because several morphological changes affected the southern slope of the cone during the activity in 2015 (Figure 1b,c). Other than this approximation, the raster operation performed between the 2012-DSM and 2017-DSM defined the lava flow extensions well (Figure 5), and the calculated volume for the south lava flow (SLF) is probably overestimated because, with respect to the 2012 topography, at least one lava flow lies below the 2016 lava flows. To calculate the error in the volume estimation (z-axis), the altitudes of control points on known morphological features that were not affected by the eruption were compared in both the 2012 and 2017 DSMs. The differences in the heights were on the order of +/−1.5 m. This value was multiplied by the area covered by each lava flow to obtain an estimate of the error of the volume calculation (Table 4). The volume of the 10–11 July 2015 BAF deposits was previously calculated [6,28], and this result is taken here to estimate the total magma volume involved in the 2014–2016 eruptive phase (Table 5).



**Table 5.** Estimated parameters for the 10–11 July 2015 pyroclastic flows emplaced along the Montegrande ravine.


**Figure 5.** Raster (7.5 m resolution) resulting from the subtraction between the 2012 and 2017 DSMs showing the main changes related to lava flow emplacement during the 2014–2016 eruptive phase. The raster was classified for values > 1 m considering the z-error of −/+ 1.5 m. The black arrow points to pyroclastic and lava products that were emplaced in 2013. (Map Projection: Universal Transverse Mercator—UTM).

**Figure 6.** (**a**) 28 November 2014 SPOT image in which the 2014 lava flow areas are clearly visible. By comparing this image with the EO-1 (ALI) (bands 3, 2, 1) of 18 July 2015 (**b**), it is possible to observe that, by November 2015, the SW lava flow was still advancing (white arrows, the dotted line points to the final extent). In this image, the lava flow emplaced just after the July eruption is also visible (black arrow). This flow was subsequently covered by the S lava flow in 2016 as evident in the 13 April SPOT image (bands: 3, 4. 2) in (**c**). (**d**) Raster showing the lava flow thicknesses as obtained by the raster operation between the 2012 and 2017 DSMs (Map Projection: Universal Transverse Mercator- UTM).

#### *4.3. Using the Etna Lava Flow Model Simulation Code*

The previously described pre- and post-eruption DSMs constitute an excellent opportunity to calibrate lava flow simulation software for the construction of hazard maps of this volcanic phenomenon.

The calibration stage is essential for the simulation of different scenarios of lava flows (or for any type of volcanic process in general) and consists of reproducing the extensions and runouts of known lava flows, which requires knowledge of the topographies prior to emplacement.

Here, we used the Etna Lava Flow Model (ELFM) code [39], which has previously been applied to other Mexican volcanoes, including Popocatépetl and Ceboruco [40,41]. The ELFM is a probabilistic model that is similar to other types of lava flow models [42,43] and simulates different paths that lava might follow when emitted from a volcanic vent [39,44]. In addition to the digital elevation model of the volcano (DEM) and the coordinates of the vent, the input parameters for the ELFM are the following: the maximum thickness of the flow (m), the possible maximum length that the lava can cross on the DEM (in the number of steps or pixels), and how the thickness of the flow varies during emplacement until its movement stops [39]. The software is based on Felpeto's algorithm [45] but incorporates improvements to solve the flow loop problem with a dynamic DEM, which makes it possible to establish the length and the height of each lava flow simulation. Moreover, these improvements allow for changes to the lava flow height over the flow path based on the distance via the use of different functions (i.e., constant, linear or logarithmic). This last parameter correlates with the viscosity and yield strength of the flow during its emplacement [46], which enables the simulation of more evolved lava flows compared with those of Etna (which are basaltic in composition) for which the software was originally created. Another variable parameter is the number of iterations in each run of the code. A value of 1000 iterations is commonly considered sufficient (i.e., no significant changes occur with greater numbers of iterations) as verified in previous studies [39,44,45,47].

In the above-mentioned research, it was necessary to reconstruct the paleo-topographies of historical lava flows because digital elevation models were not available prior to the emplacement of the flows.

This process, which is primarily performed manually, incorporates considerable error into the calibration and simulation processes because the results of simulations of lavas are highly dependent on the resolution of the DEM used [37], as repeatedly noted in simulations of other volcanic and geological flow phenomena [48–51]. For the case of Volcán de Colima, the September 2014 and September 2016 lava flows were reproduced on the existing DSM obtained from the 2012 INEGI LIDAR dataset [37], which was resampled to 10-m horizontal resolution.

The results of the ELFM were ASCII (see glossary) files that were processed in the ArcGIS 10.2® software according to the following steps: (a) transformation the ASCII files to raster files; (b) reclassification of each raster file into eight groups of values, excluding the pixels with values between 1 and 5 (which corresponded to 0.5% of the total iterations and were considered to be overestimations in the simulations [44]); and (c) transformation of the raster files into shape files (polygons) to calculate the inundated areas in each of the calibration runs and compare them with the real areas of the considered lava flows.

#### **5. Results**

#### *5.1. Lava Flow Inundation Area and Volume Estimation*

Three main lava flows were identified from the analysis of the SPOT images, and they are named here based on the directions of their emplacements as follows: May 2014 West Lava Flow (WLF), September 2014 South-West Lava Flow (SWLF), and September 2016 South Lava Flow (SLF) (Figure 6 and Table 4). Based on the raster operations, their volumes were estimated and ranged from a minimum of 7.5 <sup>×</sup> <sup>10</sup><sup>6</sup> m3 for the WLF to a maximum of 16 <sup>×</sup> <sup>10</sup><sup>6</sup> m3 for the SLF (Table 4). As previously mentioned, the volume estimation for the most recent SLF is probably overestimated because the 2012-DSM does not reflect all of the morphological modifications of the southern slope of the cone during and after the 10–11 July 2015 eruptive activity. Lava flows were emplaced over slopes of 29.5◦ for the WLF, 27.8◦ for the SWLF and 31◦ for the SLF and reached maximum distances of 1.2, 2.2 and 2.2 km, respectively (Table 4). Interestingly, in the SPOT images from 28 November 2014, the SWLF is still advancing because its extent does not correspond with the maximum runout visible in the 18 July 2015 EO-1 (ALI) and 2017 SPOT images (Figure 6a–c). Therefore, by 28 November 2014, the lava traveled 1.9 km (for a volume of 4.7 <sup>×</sup> 106 m<sup>3</sup> based on the 2012- and 2014-DSM raster operations; Table 4). Its final runout is 2.2 km for a total volume of 12.2 <sup>×</sup> 10<sup>6</sup> m3 (Figure 6a,b). For this lava flow, the exact dates of the initial and final extrusion times are unknown. If its initial stage of extrusion was detected in September 2014, by 28 November (~75 days), an extrusion rate of ~0.7 m3/s could be estimated, which is in the range of the effusion rate calculated by [52] for the same lava flow. If this estimation is correct, the lava flow probably stopped by March 2015, but no images were found to confirm this assumption. In contrast, a more precise estimation could be performed for the SLF. A thermal anomaly related to the SLF was first reported on 27 September 2016, and the lava flow maximum runout was reported on October 19th [23,24]. Based on this information, an extrusion rate of ~8 m3/s was estimated here, which is much higher than that of the SWLF but very similar to the extrusion rates calculated for lava flows in 2004 [53]. From the raster values of the lava flow bodies and based on the topographic profiles, the lava flows exhibit a maximum thickness of approximately 50 m towards their distal fronts (Figure 6d). The WLF shows a general convex profile with a clear blocky superficial texture, which contrasts with the SWLF, which shows a concave profile with a lateral levee (Figure 6; Figure 7). The SLV consists of three main lava bodies with concave profiles similar to the SWLF.

**Figure 7.** Longitudinal and transverse profiles of the WLF, SWLF and SLF described in the text. The transverse profiles were delineated in the middle part and at the distal end of the path of each lava flow. (Map Projection: Universal Transverse Mercator—UTM).

#### *5.2. Block-and-Ash Flow Deposits: Inundation Limits, Volume and Thickness Variation*

The 10–11 July 2015 collapse of the summit dome produced pyroclastic flows that inundated the Montegrande ravine (Figure 2). The 2015 SPOT product was segmented to define the affected area [6] (Figure 2b) and distinguish between the valley-confined facies that consist of massive, matrix-supported BAF deposits that are up to a 20 m thick with cm to m blocks imbedded in an ashy matrix, and the overbank facies that correspond to 1 to 3 m thick massive BAF deposits that are covered by centimetric layers of massive ash. These fine, homogenous ash layers on top of the overbank unit comprise the feature that allowed for its discrimination from the valley-confined facies during image segmentation. A volume of 4.2 <sup>×</sup> 106 m3 was estimated for the valley-confined BAF based on the difference between the 2012-DSM and the 2015-DSM (Table 5). The dense rock equivalent (DRE) [54] was then estimated based on an andesitic lava density of 2.4 g/cm<sup>3</sup> and a BAF deposit density of 2 g/cm3, yielding a DRE (magma) volume of approximately 3.5 <sup>×</sup> 106 m3. [28] analyzed two 50-cm resolution stereopairs of Pleiade-1A satellite images acquired in April 2013 and January 2016. Based on photogrammetry of these two stereo pairs, these authors obtained two 1 m spatial resolution DEMs based on which a volume of 5.8 <sup>×</sup> 10<sup>6</sup> m3 for the same valley-confined facies was estimated. Considering the errors (Table 5), both volume estimations are comparable.

#### *5.3. Calibration of the ELFM at the Colima Volcano*

Software calibration is a fundamental step prior to lava flow simulations of different eruptive scenarios for hazard evaluation. It consists of several test-error simulations and modification of the input parameters until the best fit of the area covered by the real lava flow is found. The 2014 SWLF and the 2016 SLF were reproduced on the 2012 DSM [37] after totals of 42 and 25 runs, respectively, using 1000 iterations in each run as a fixed parameter.

A fundamental aspect of the calibration was the selection of the exact possible vent of each of the lava flows. Tests were performed on eight different pairs of coordinates for the 2014 SWLF and three different vent position for the 2016 SLF, all of which were located in the proximal area of the lava flow polygons obtained by the raster subtraction (Figures 5 and 6d). Notably, the positions of the three emission centers used for each flow are separated by an average of 30 m.

For both lava flows and each of the vents selected, tests with linear and logarithmic variations of the thickness of the flow were performed and found that a linear increase in the thickness (value 1 in the software) over the path of the flows was the best choice for the lavas emplaced on high slopes; this is also the case for the 2014 and 2016 Colima lava flows (slopes between 28◦ and 31◦, Table 2) in which the maximum thicknesses (35 to 50 m) are reached in the middle and distal portions of the total lengths (Figure 7).

The other input parameters needed for the ELFM were varied in each of the runs to search for the datasets that best fit the extents of the real lava flows. The maximum thickness of the flow was modified between 10 and 40 m, and the possible maximum length (in number of steps or pixels) that the lava could cross was varied between 300 and 550 steps. The best parameters for reproducing the Colima 2014 SWLF were a thickness of 30 m and 400 steps (Figure 8a). For the 2016 SLF, a thickness of 30 m and 370 steps were sufficient to reach the length of the observed lava flow (Figure 8b). A thickness value of 30 m was selected as the average value of the thickness of the lava flows. The maximum value of 50 m, which was obtained in the calculation of the thickness in the raster operations (Figure 4d), was not considered in the calibration simulations because this value was only obtained at specific points of the total area of each lava flow, as observed in the transverse profiles (Figure 7).

The selected simulations revealed a low percentage of underestimated areas (i.e., areas flooded by the real lava flows and not covered by the simulation) and fit with the full spatial extent of the reproduced lavas. However, in both cases, the simulations overestimated the flooded areas—by 40% in the case of the SWLF and 35% in the case of the SLF. This overestimation was probably due to the fact that the 2012 DSM does not account for changes at the summit of the cone that occurred between 2012 and 2014.

**Figure 8.** Results of the simulations showing the best fits with regard to the extensions of the real lava flows. The simulations were performed with 1000 iterations. The ASCII (see glossary) files resulting from the Etna Lava Flow Model (ELFM) software runs were transformed into raster files and later reclassified in a discrete number of groups of iterations as shown in the legend. (**a**) Reclassified raster of the SWLF simulation using 400 steps for the flow path, 30 m for the final thickness and a linear increase in the flow thickness. (**b**) Reclassified raster of the SLF simulation using 370 steps for the flow path, 30 m for the final thickness and a linear increase in thickness. (Map Projection: Universal Transverse Mercator- UTM).

Nevertheless, based on this calibration, a lava flow hazard evaluation could be extended to the entire cone, including other vents around the crater rim, for lavas with similar extensions and volumes, and also for other eruptive scenarios that could involve the emission of lava flows of greater length and volume.

#### **6. Discussion**

Using SPOT and EO-1 (ALI) data, we were able to map and estimate the volume of the deposits emitted during the 2014–2016 eruptive phase of Volcán de Colima, which is one of the most active volcanoes in Mexico. These data were freely available, and their temporality acquisition was optimal for defining the chronology of the eruption and the maximum extensions of the pyroclastic flows and lava flows, their volumes and their extrusion rates. Although products of higher resolution are available, such as the Pleiades 1A data (Pleiades-HR is composed by two-spacecraft constellation, optical high-resolution panchromatic -0.7 m- and multispectral -2.8 m-), most are prohibitively expensive, especially for institutions working in Latin America. Despite this limitation, the results presented here demonstrate that SPOT images can be a useful tool and provide a spatial resolution that is sufficient to obtain DSMs from which volume estimations can be retrieved. Moreover, these volume estimations agree with those obtained, for example, with DEMs based on Pleiades 1A data [28] for the BAF deposits (Table 5) and by direct observation of lava flows [52]. A map of the total affected area that distinguishes between the lava flows and pyroclastic deposits is presented here (Figure 9). This map is an invaluable tool that can be used to calibrate numerical simulations and upgrade the hazard map of Volcán de Colima, for which a similar scenario has not previously been considered [9]. Specifically, the calibration of the lava flow simulation software presented here benefited from the availability of the real topography prior to the event and the total distribution of the simulated lava flow. Based on the best-fit values obtained in this work, a range of values for each parameter can be established to

simulate the emission of andesitic lava flows of a few kilometers in length from eruptive centers on the edges of central crater and to establish other sets of data to reproduce different scenarios of more voluminous lavas.

**Figure 9.** Map showing the distribution of lava flow deposits and PDC (Pyroclastic Density Currents) deposits associated with the 2014–2016 eruptive phase of Volcán de Colima. (Map Projection: Universal Transverse Mercator—UTM).

An additional result of this study relates to the eruptive behavior at Volcán de Colima. Based on the chronology of the effusive phases and the volume of the lava flows, the extrusion rates were estimated. The calculated extrusion rates spanned from a minimum of 0.7 to a maximum of 8 m3/s, in agreement with the information reported in a previous study (i.e., [52]). The lower rate corresponds to the SWLF that emplaced over a lower slope gradient (27.8◦) and with a lower magma volume (12.2 <sup>×</sup> <sup>10</sup><sup>6</sup> <sup>m</sup>3), and the higher value corresponds to the SLV that emplaced over a 31◦ slope with a larger volume (16 <sup>×</sup> 10<sup>6</sup> m3). The correlations between the slope gradient, magma volume and extrusion rates fit with the general behavior of lava flow [55], which again confirms the validity of the parameters calculated for the 2014–2016 lava flows.

Finally, based on the estimation of the total volume of the magma involved in the 2014–2016 eruptive phase, which was considered to be an anomalous scenario associated with historical activity (i.e., [6,7]), some considerations about the eruptive style can be presented. Considering the lava flows and pyroclastic flow deposits, we estimated a total magma volume of approximately <sup>40</sup> <sup>×</sup> <sup>10</sup><sup>6</sup> m3. This value is much higher than previous activities associated with dome collapses, which involve BAF volumes on the order of 10<sup>5</sup> m3 (i.e., [4,5]), but one order of magnitude lower than that of the 1913 Plinian eruption, which had a total volume of 3 <sup>×</sup> <sup>10</sup><sup>8</sup> m3 [8]. Thus, even if the 2015 scenario was unexpected (primarily due to the maximum runout of the BAFs, which have usually extended up to 7 km during previous eruptive episodes), the magma volume in the magma chamber was at the lower limit of the estimation at which a Plinian or even sub-Plinian scenario could be expected based on numerical models [56]. A larger magma volume (108 m3) is needed to generate an overpressure sufficient to trigger a Plinian scenario. This is an important result that helps to better understand the eruptive behavior of Volcán the Colima and supports the theoretical results obtained with numerical models.

#### **7. Conclusions**

This analysis focused on defining the chronology of the events, the affected area and the total magma volume involved in the eruption. As demonstrated here, the SPOT data were optimal for generating a first approximation of the magnitude of the eruption and producing a map with the distribution of the deposits. These are invaluable data for understanding a volcano's behavior and simulating possible future scenarios, especially considering the reduced accessibility of the majority of the area around an active volcano. Based on here presented results, the hazard assessments for the block-and-ash flows and lava flows at Volcán de Colima can be revised and improved. Additionally, based on numerical models [56], the estimated magma volume involved in the 2014–2016 eruptive phase was too low to trigger a Plinian scenario similar to the 1913 activity. A larger magma volume stored in the magma chamber would be needed, and, as observed in the precursory activity of the 1913 eruption, dome destruction would be accompanied by explosive events; no such events were observed during the 2015 activity.

Remote sensing is an invaluable tool for the rapid mapping of active volcanoes after major eruptions. As recently observed at Volcán de Fuego in 2018 (Guatemala), a quick mapping to evaluate damage can be realized and used for rescue activities in the few days following a disaster (i.e., Information Technology for Humanitarian Assistance, Cooperation and Action (ITHACA), http://www.ithacaweb.org/; United Nations Institute for Training and Research Operational Satellite Applications Programme (UNITAR-UNOSAT), www.unitar.org/unosat [57]). Subsequently, the same products can be used to define a map of the different deposits (lava flows and pyroclastic flow) that enable volume estimation, identify key features for numerical model calibration, and allow the construction of hazard maps.

**Author Contributions:** Conceptualization, L.C. and D.F.; Data curation, N.D.; Formal analysis, L.C.; Investigation, J.C.G.-R.; Methodology, N.D., L.C., D.F. and P.F.; Software, P.F.; Supervision, L.C.; Writing—review & editing, N.D., L.C., D.F. and J.C.G.-R.

**Funding:** This research was funded by the PAPIIT-DGAPA project IN105116 to Lucia Capra. DF was benefited by a postdoctoral fellowship from the Conacyt project PN-360 to LC.

**Acknowledgments:** We thank to Station of Telemetry Reception in Mexico (Estación de Recepción México-ERMEX) for having provided all stack of SPOT6/7 images. DF thanks Dr. Gianluca Groppelli, of the Istituto per la Dinamica dei Processi Ambientali (Milan, Italy), who provided access to ELFM software and advised the work of modeling lava flows. Finally, we thank to earthexplorer webserver (https://earthexplorer.usgs.gov/) for having provided EO-1 (ALI) image.

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