1. Introduction
Explosive volcanic eruptions commonly produce one or more products (e.g., tephra fallout, pyroclastic flows, lava flows/domes, etc.) that pose risks to the lives and livelihood of local populations. These products can threaten thousands of lives in highly populated regions [
1]. For example, pyroclastic density currents (PDCs) have caused ~90,000 fatalities historically, the most of any volcanic hazard [
2]. PDCs are high-temperature (200–800 °C), high-velocity (80–750 km/h) mixtures of rocks, ash and gas that typically have runout distances between 10–15 km, but in some cases, as long as 100s of kilometers [
3,
4]. PDCs are comprised of a dense and dilute component. The latter is not confined by topography and is known as either a pyroclastic surge (PS) or a dilute PDC [
5,
6]. The dense portion, also known as a pyroclastic flow (PF) or dense PDC, tends to follow existing topographic lows. The formation of PDCs commonly occurs as the result of either the collapse of a large convecting eruptive column, a “boiling-over” of PDC material over the crater rim, or the presence of a lava dome that decompresses from gravitational collapse, or serves to increase the magnitude of an explosive eruption following its emplacement.
Because of the extreme hazard potential of PDCs it is critical to investigate volcanoes that have a history of producing these events, one recent example is Fuego volcano in Guatemala. The June 2018 eruption of Fuego produced a large PDC, resulting in the mass casualty event. Approximately 50,000 people live within 10 km and ~1,000,000 within 30 km of the summit of Fuego [
7]. Numerous communities (e.g., El Rodeo, La Reina, Los Lotes) are located within ~5 km and are at the highest risk. However, any individuals living within ~8 km are at risk from a PDC due to the long runout distances at this volcano because of the volcanic geomorphology [
8]. Historically, PDCs and other eruptive products have been confined to Fuego’s large drainage ravines (barrancas) that are distributed radially around the summit in all but the northerly direction, as shown in
Figure 1 [
9]. However, if these flows leave the confines of the barrancas (as happened in 2018), the threat increases significantly to a larger area and more people.
Over the past five decades, the quantity of data available on volcanic eruptions via satellite-based remote sensing has increased exponentially [
10]. Multiple studies have made use of these data to investigate specific events or styles of activity (e.g., [
11,
12,
13]). These studies have either examined volcanic activity globally or focused on a specific geographic region or individual volcano. Several more recent studies (e.g., [
7,
14,
15]) have used remote sensing to investigate the eruptions, emissions, and eruptive history at Fuego, specifically, utilizing a variety of satellite remote sensing instruments, such as the Moderate Resolution Imaging Spectroradiometer (MODIS), the Infrared Atmospheric Sounding Interferometer (IASI), and CubeSats. These have helped to fill a knowledge gap that existed for this volcano; specifically, the increase of volcanic activity since 1999.
Our study uses data from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) instrument. ASTER is one of the five instruments on the National Aeronautics and Space Administration (NASA) Terra satellite [
16]. These data have been used to monitor and assess volcanic activity during numerous past eruptions (e.g., [
11,
12,
17]). ASTER has a 15 m per pixel spatial resolution in the visible/near infrared (VNIR), which is also used to create digital elevation models (DEMs), and a 90 m resolution in the thermal infrared (TIR). The TIR sensor has a noise equivalent delta temperature of ~0.15 °C, allowing derivation of accurate temperatures following atmospheric correction and temperature/emissivity separation of the radiance-at-sensor data (e.g., [
18]). These data saturate at pixel-integrated brightness temperatures in excess of 100 °C.
The ASTER data are critical to assess the directional tendency of PDCs over the past 15 years. The volcanic process and product causing each thermal anomaly detected by ASTER were determined and then validated by comparing them to ground-based reports. This allowed the number and direction of all PDC events to be determined. The ASTER DEM data also provided the pre-flow topography for the subsequent modeling using VolcFlow. The ultimate goal of this study was not to refine a flow model that perfectly matched the 2018 deposit, but rather to assess how such numerical flow model freeware performs using free, globally-available DEMs of varying quality. This provides a basis from which future studies at other volcanoes could be assessed. We performed modeling at Fuego on two different ASTER-DEM products to quantify the DEM quality errors, as well as the most common direction of future events. The historical TIR data analysis, therefore, directly informed the flow modeling of future PDCs. The combination of remote sensing with a DEM uncertainty analysis is unique and important for the creation of a PDC-specific hazard map, which shows that the greatest threat is to the populations nearest the summit and those closest to the barrancas.
1.1. 3 June 2018, Eruption
The 3 June 2018, eruption of Fuego produced a large explosive phase resulting in a column collapse and generation of a PDC, which also incorporated material from a partial collapse of the SE flank. Unlike previous PDCs that commonly remain confined to the barrancas, the volume of this PDC caused it to leave the confines of the Las Lajas barranca in two locations (white stars in
Figure 1). The first at ~5 km from the summit, and the second ~3.5 km further downslope. The first exit point devastated the La Réunion Resort (
Figure 2). There were no reported fatalities at this location, because of the successful evacuation of the resort staff and guests prior to the PDC emplacement. The second exit point was ~1.5 km to the north of the village of Los Lotes (
Figure 3). Unlike the La Réunion Resort, Los Lotes was not evacuated and experienced the majority of casualties.
Official incident reports issued by the United Nations Office for the Coordination of Humanitarian Affairs (OCHA) determined 113 were killed from the event (
Table 1). News reports following the incident placed the death toll closer to 200. As of the final OCHA report covering 19–27 June 2018, there were still 197 people missing. However, the most recent study documents more than 400 deaths [
19].
1.2. Background
Fuego volcano is located in central Guatemala and has a well-defined summit crater at an elevation of 3763 m. It forms half of the Fuego-Acatenango volcanic complex and is located within the second of eight segments of the Central American Volcanic Arc [
20,
21,
22]. Fuego has experienced five VEI = 4 eruptions in 1581–2, 1717, 1880, 1932, and 1974 [
8,
23,
24]. During its history Fuego has also had more than 60 subplinian eruptions and longer periods (months to years) of low-level strombolian volcanism [
20]. This more common activity occurs many times daily and has little to no effect on the local people living in close proximity to the volcano. Since 1999, a new eruptive style began that consists of persistent strombolian activity, lava fountaining, and explosions [
20]. It is punctuated by the occasional paroxysmal eruption of greater energy that can produce lava flows and PDCs [
8]. Naismith et al. (2019) proposed that starting in 2015 Fuego entered into a different period of even greater volcanic activity. This period follows eruptive activity described by Lyons et al. (2009) but is characterized by a greater frequency of paroxysmal eruptions with the following series of events: (i) effusion of lava flows and an increase in summit activity, (ii) an intense paroxysmal eruption lasting 24–48 h that produces a sustained eruptive column and PDCs, followed by (iii) a decrease in activity [
7]. Our study focuses on the PDCs that formed during one of these paroxysmal eruptions on 3 June 2018.
3. Results
For the study period of 31 January 2003, to 30 April 2018, 193 nighttime ASTER TIR scenes were found. Of those, ~57% (109) had thermally-elevated regions above the scene’s background temperature, whereas the others either had cloud cover that obscured the summit or had no detectable activity. Within the subset that had detected activity, 47 were single to multipixel summit-confined hot spots, and 62 showed evidence of a flow.
Figure 4 shows examples of four volcanic events: a lava flow, a summit thermal anomaly, and two positively identified PDCs confirmed using INSIVUMEH reports.
After the analysis of each scene and comparison with the INSIVUMEH reports, a total of 54 PDCs were confirmed during the ~15-year study period. The average flow length measured from the ASTER data was ~5 km. Their directional distribution is dominantly southward as might be expected based on the topography of the Fuego-Acatenango complex (
Figure 5). The largest number of events (16) traveled southwest down the Taniluya and Cenizas barrancas. Due to the spatial resolution of the ASTER TIR data, it is not always possible to differentiate between these two barrancas; in which case, the INSIVUMEH reports were used for validation. The remaining PDCs were nearly equally distributed to the northwest, south, and southeast, with very little activity in the Honda barranca to the northeast.
The VolcFlow model runs correlate with the historical reports and findings of previous studies, which indicate PDC events generally remain inside the barrancas as they were emplaced [
9]. This is expected because PDCs at Fuego are generally more PF-dominated as opposed to PS-dominated. The observation that most PDC’s from Fuego are dominated by dense material is also supported by a lack of PS deposits seen in the remote sensing data or reported by INSIVUMEH.
VolcFlow was first run using the 30 × 10
6 m
3 volume estimate of the 3 June 2018, PDC and using the ASTER GDEM product for the topographic base. A starting location/elevation was chosen to coincide with the head of the Las Lajas barranca, rather than directly at the summit. The PDC flow coverage for this initial result is shorter and distributed in more barrancas than what was observed following the June 2018 eruption (
Figure 6). Furthermore, the model was unable to replicate the PDC leaving the barranca at the two locations during that event.
A single-scene DEM was chosen to assess how the higher amounts of pixel to pixel noise (compared to the GDEM) would affect the model results. By the ASTER day overpass time (~10:45 am local time), higher percentages of clouds are commonly encountered over the mountainous regions of Guatemala, even during the winter dry season. The formation of these orographic clouds typically begins about an hour earlier. Therefore, it is difficult to find scenes that are cloud-free. The scene closest to the eruption date with the least cloud cover was acquired on 24 January 2018 (
Figure 6b). The DEM generated from this scene has the typical minor topographic errors caused by the pixel to pixel noise, as well as dense clouds north of the summit and much smaller clouds elsewhere (
Figure 7). The presence of dense clouds to the north create a large false topographic high over Acatenango, which did not affect the VolcFlow results. The pixel noise, however, did result in modeled flow paths and lengths that were inconsistent compared to the GDEM model run. Use of the DEM with VolcFlow resulted in ~28% less area covered by the deposit and flow lengths that were ~54% shorter than the results using GDEM.
To quantify the impact of this DEM topographic error on the model, the single-scene DEM data were subtracted from the GDEM data (
Figure 7). The minor cloud cover around the summit created positive topographic errors as high as 475 m above the baseline, whereas shadows from these clouds caused negative errors as low as 238 m below. To determine the average difference between the GDEM and the single-scene DEM, four 1 km
2 regions of interest (ROI) were chosen in areas of relatively consistent terrain that were free of clouds and cloud shadows. The results are shown in
Table 3.
The pixel/scan line noise present in the single-scene DEM did cause artificial elevation changes that resulted in modeled flow paths that were shorter and with different routings. Although the average elevation change in the ROIs is within the estimated error, the GDEM contains less random noise and more clearly defines the barrancas. The minor amount of cloud cover in the single-scene DEM, however, did not have a significant impact on the individual model run. This despite the clouds/shadows producing the extremes in the difference image, which are significantly outside ASTER DEM error estimates [
30].
VolcFlow was next run using the GDEM for three increasing volumes to produce a generalized hazard map (
Figure 8). For this facet of the analysis, the summit location and elevation were used for the model starting point. The hazard map tends to follow topographic lows as expected. Slight deviations from these occur where the forward momentum of the modeled flow approaches zero and/or where the flows extend beyond the confines of a topographic low. Model results with higher eruptive volumes do increase the potential risk to communities further from the summit of Fuego, as would be expected. The maximum modeled runout distance for a 40 × 10
6 m
3 volume flow is ~10 km from the summit, for example. However, even the smallest volume produces modeled flows that reach the two population centers (Los Lotes and the La Réunion resort) that were directly impacted by the 3 June 2018, eruption, confirming the usefulness of this simplified approach (
Figure 8).
4. Discussion
Satellite-based remote sensing is an effective method for observing synoptic volcanic activity, including eruptions and eruptive precursors (e.g., [
12,
17,
41]). However, sensors with the higher spatial resolution needed to resolve smaller-scale details commonly have poorer temporal resolution, which can result in events being missed. In contrast, those with a higher temporal cadence, such as MODIS or the Advanced Very-High-Resolution Radiometer (AVHRR), needed for capturing all phases of an eruption, have much lower spatial resolution that precludes resolving details, such as accurate flow dimensions and direction. Direct ground-based observations (where possible) of volcanic deposits following an eruption do provide a much greater understanding of the event (e.g., [
42,
43]). However, where continuous ground observations are not feasible, reliance on local monitoring groups/agencies and activity reports are helpful. Having ground reports available from INSIVUMEH to validate the remote sensing data used here, for example, was important for this study. However, both the ASTER data and INSIVUMEH reports do have limitations.
ASTER has an average temporal resolution of ~8 days for Fuego, which is higher than the reported value of ~16 days [
16]. The higher volume of ASTER TIR data was possible because of the presence of the ASTER Urgent Request Protocol (URP) program [
44,
45]. The URP integrates ASTER into a sensor web architecture, being triggered by thermal detections from MODIS. These triggers, which happen continuously during a larger eruption, produce dozens more scheduling requests for ASTER data. By providing higher temporal, high spatial resolution TIR data, the URP program becomes an important tool for monitoring volcanic trends. However, even with this improved temporal cadence, volcanic events can still be missed entirely and/or only partially captured if an eruptive product (e.g., lava flow or PDC deposit) cools to become undetectable in the TIR wavelength region. We believe this influenced our PDC measured lengths somewhat resulting in slightly shorter average values.
The ability to remotely observe these volcanic products in the TIR data at a high spatial resolution is critical to this study. Without the high spatial resolution and accurate radiometric resolution of the ASTER TIR data, it would be impossible to determine the flow direction and deposit type produced from each eruption. The only other publicly available TIR data near this spatial scale and acquired over a long time period are the data from the Landsat satellites. In the past, the Landsat TIR data were not routinely acquired at night, although this is becoming more common for volcanic and fire activity. The Landsat TIR data do not have the ability to point off-nadir as ASTER TIR does, which decreases the time between observations. Furthermore, the limit of one or two TIR spectral bands, depending on the Landsat TIR instrument design and/or the health of that instrument, fundamentally reduces the accuracy of derived temperatures from those data. Despite those limitations, a long-term TIR study, such as this one, conducted at other volcanoes, could be augmented by Landsat data over the same time period studied or extended earlier in time using Landsat TIR data acquired before 2000.
The Naismith et al. (2019) study concluded that starting in 2015 Fuego entered a new period of increased activity that produces more lava flows and PDCs [
7]. This was partially concluded using the MIROVA monitoring system, an automatic volcanic hotspot detection system based on the 1 km TIR data from MODIS [
46]. Based on our analysis of the ASTER archive, we see a ~50% rise in the number of scenes containing elevated thermal data (e.g., lava flows, PDCs, and lahars) beginning in 2015, which supports the findings of Naismith et al. (2019).
Satellite remote sensing of volcanic activity is not without limitations and whenever possible should be collaborated with ground reports. However, INSIVUMEH reports also have limitations. Nine events identified as potential lava flows or PDCs were observed with ASTER but had no corresponding INSIVUMEH report. This could be due to a number of factors, such as the summit being obscured by clouds or the event being too small for ground observers to see. Furthermore, on the morning of the 2018 eruption, authorities announced that the volcano “looked fine”. These human-induced limitations emphasize the usefulness of data synthesis between satellite and human monitoring at Fuego, each at different time and spatial scales, to provide a much more accurate chronology of the volcanic activity over the study period. The nearly 15 years of ASTER TIR data of Fuego volcano compared to and validated by the INSIVUMEH ground reports indicate that the most common PDC flow direction is to the south, closely followed by flows to the southwest and southeast. This information, combined with continued TIR data and DEM collections, can be used in future developments of disaster response and mitigation plans for the local monitoring agencies.
VolcFlow can model PDCs, avalanches, and lava flows assuming the material property values (e.g., density and cohesion) are available and accurate. These properties are difficult to calculate directly from remote sensing data. For this study, values from prior VolcFlow modeling studies of PDCs were used; however, these could be refined with laboratory analysis of samples from the June 3 flow. The estimated volume of the PDC used in the study came from Naismith et al. (2019) and is consistent with previous eruptions of Fuego [
8]. However, a thorough field campaign and more accurate pre-eruption DEM would provide more accurate constraints on the volume of material erupted.
A second important model input for VolcFlow is the pre-flow DEM [
35]. For this study, two versions of DEMs generated by the ASTER instrument were used. The stated vertical accuracy of the ASTER GDEM product is 20 m, better than a single-scene DEM accuracy of 30 m. The GDEM achieves this improved accuracy and reduction in pixel to pixel noise through averaging large numbers of single-scene DEMs acquired over time. This can create constraints for flow modeling. First, version 2 of the GDEM uses scenes from 2000–2011. Therefore, any topographic changes within this time period are lost. However, if no major topographic changes have occurred, as with Fuego, its use is appropriate. A version 3 GDEM was recently released which uses data from 2000–2013. It was released too late to be incorporated into this study, however the differences between it and version 2 at Fuego were assessed to be negligible. Version 3 contains only minor fixes, such as filling of holes. Second, to account for any significant topographic change, a model requires the most recent (and most accurate) DEM. The use of a single-scene DEM closest in time to an event is potentially more appropriate than using the GDEM. However, as seen by our analysis of the single-scene DEM, the amount of noise and topographic anomalies can adversely impact the modeling in certain locations (
Figure 6b).
The ultimate goal of this study was not to create a model output that precisely matched the June 2018 deposit, nor assess whether VolcFlow is the appropriate flow model for this particular PDC emplacement. To achieve a more accurate fit with the VolcFlow model, its input parameters could be varied iteratively, and the subsequent output assessed. Rather, we wanted to examine the accuracy of the different ASTER-derived DEMs and their appropriateness in future modeling efforts. For this purpose, the VolcFlow model was chosen and run with a set of standard input parameters used in prior studies and the results compared using those two DEM products.
VolcFlow modeling using the GDEM was able to replicate the flow length, however there were notable differences. These are likely due to the constraints of the model to fully replicate all components of a PDC, the chosen starting location for the modeling, as well as the resolution of the DEM itself. For example, the model results show the deposit being emplaced into several barrancas, yet nearly all the material was concentrated into one barranca during the 2018 eruption. This higher volume emplaced into one barranca is likely the reason that the PDC overran the barranca walls in two locations downslope. Progressive infilling of the barranca by PDC material changed the depth in near real time, which was not resolved in the DEM. In addition, initiating the model at the summit (with a starting velocity of 0) is likely not accurate if the PDC initiated at a higher elevation from partial column collapse. The modeled results are therefore shorter (as we see in our results). Model results using the single-scene DEM presented other problems. Noise variations in the data resulted in minor random elevation anomalies from one pixel to the next that directly impacted the flow’s progression in the modeling. This produced results showing the PDC material dispersed laterally over a larger area higher on the slopes before becoming confined to multiple barrancas. The PDCs, therefore, did not travel as far as the population centers affected by the 3 June 2018, eruption. Although not a factor for this study, more extreme topography errors caused by the presence of cloud and cloud shadow in a scene would also affect the model output even more and should be considered where using individual ASTER scene DEMs for future studies. If, for example, there has been significant topographic changes at the volcano under study, the GDEM would not be an appropriate choice for a flow model. Although not employed here, some of the random noise present in these DEMs can be lessened with image smoothing and noise removal algorithms making the single scene DEM potentially more useful.
Volcano-wide hazard modeling of PDC-forming eruptions shows that population centers in close proximity to the barrancas to the southeast of Fuego (e.g., Los Lotes, San Jose, and El Rodeo) are at higher risk during larger paroxysmal eruptions. During the 3 June 2018, eruption, a PDC formed, flowed down Las Lajas barranca, and later exited it at two locations destroying both the La Réunion resort and the community of Los Lotes. These two locations are within the modeled flow hazard map for an eruption similar in size to the 3 June 2018, event. Because of the barranca-exiting risk of future large PDCs, any population center within 2 km of a barranca does have a higher risk during eruptions. However, it is outside the current ability of the VolcFlow model to recreate these more complex barranca-departure events. The most likely scenario is that these barrancas become infilled down-valley either by ongoing erosion during the rainy seasons or from the deposition of new material during the early stages of an eruption, thus changing the depth profile. Modeling future hazards at Fuego will benefit from more accurate and more frequently updated DEMs. For example, dedicated field campaigns following eruptions and flow emplacement using either drone-based LIDAR or structure from motion images could create updated DEMs over the areas of change, which could be used to augment the larger ASTER GDEM for more accurate modeling of future flow hazards.
5. Conclusions
Orbital ASTER TIR data refined using summary reports from INSIVUMEH provided detailed information on the recent PDC travel directions at Fuego volcano. Although the reports provided insights into the deposit type and specific barranca occupied, the resolution of the ASTER data over 15 years, resolved the activity and specific flow events. This study would have been less accurate without either the reports or the high spatial resolution TIR data acquired with improved temporal frequency enabled by the ASTER URP Program.
Despite the June 2018 PDC-forming eruption at Fuego occurring seven years after the data used to create the GDEM, it was found to be more appropriate and accurate for modeling the general trends of the flow produced by that eruption. Any topographic changes that did occur from 2011 to 2018 were small enough in comparison to the increased noise present in the single-scene DEM. However, this may not be the case at other volcanoes, especially those with significant topographic changes since 2011. A similar approach is easily adaptable to other volcanoes where the observer reporting may or may not be available and/or incorporating TIR data from Landsat. For any modeling effort integrated into those studies, the choice of flow model and the DEM are important factors. The ability of ASTER to provide several DEM products globally and over time make it the best option for most cases. However, the differences in the DEM products (e.g., the single-scene versus the GDEM) and the impact on the flow model’s uncertainty must be considered.
Despite the lack of a precise match between the model output and the June 2018 deposit, the historical tendency of flow directions as seen in the ASTER data analysis supports the VolcFlow model results, which identifies the risk to the surrounding population centers. The conditions leading to the two barranca-departure events were beyond the modeling and the DEM resolution. These events are what caused a majority of the fatalities and therefore require further analysis with more sophisticated modeling and progressively updated, higher-resolution DEMs. Despite this limitation, the generalized PDC hazard map with a PDC volume similar to that of the June 2018 volume identifies the population centers of La Réunion and Los Lotes as at-risk areas. With Fuego possibly now in a period of increased volcanic activity, the results of this study are important to potentially help mitigate future PDC disasters. As ASTER ages, similarly-conducted studies will benefit from a continuation of future TIR data with specifications equal to or improved from that of ASTER.