An automated ash dispersion forecast system: case study Popocatépetl volcano, Mexico

Abstract

An operational volcanic ash dispersion forecast system was developed for Popocatépetl. It runs automatically every day developing 108 possible scenarios of ash dispersion for the following 36 h. Scenarios are simulated for three eruption column heights: 3 km, 5 km, and 10 km above the volcano’s crater level, every hour for eruptions lasting 1 h. For each hypothetical eruption that starts every hour, the dispersion during the following 8 h is modelled. The system uses the Weather Research and Forecasting (WRF) model for weather data and the Fall3D model. It includes a visualization website that displays, among other products: ground accumulation, deposit load, and concentration at relevant flight levels. Popocatépetl volcano, located ~ 60 km from Mexico Megacity was selected as a case study. A comparison from ash forecast system results and satellite observations is presented. The system developed and tested here can be adapted to be operative at any volcano.

Introduction

Ash transport and deposition can affect economic activities (e.g. civil aviation routes, airports and land routes), infrastructure (e.g. roof collapse due to ash load), cause health effects (e.g. respiratory problems, eye and skin irritation) and other environmental impacts (water contamination, animal deaths, agricultural and forestry impacts, etc.) (Montiel et al. 2022; Bia et al. 2022; Elissondo et al. 2016; Rivera-Tapia et al. 2006; Horwell and Baxter 2006).Volcanic ash can take from minutes to hours before reaching and deposit over dense urban areas (Macias et al. 1995). The impact area from volcanic ash emissions are issued by the Volcanic Ash Advisory Centres (VAACs) by means of satellite imagery and simulations from Volcanic Ash Transport and Dispersion (VATD) models (Harvey et al. 2020).

Efforts around the world for the forecasting of volcanic ash clouds have been evolving to prevent hazards. For instance, at Etna volcano, Italy, a system was developed based on monitoring (e.g., remote sensing) data to forecast ash dispersal and accumulation, and then to produce and publish hazard maps to be consulted by civil protection authorities (Scollo et al. 2009). In such monitoring/forecasting systems, the flow of information has a considerable delay from the eruption event to the ash dispersion forecast. Rizza et al. (2020) used the Weather Research and Forecasting model coupled with Chemistry (WRF-Chem) for the simulation of volcanic ash and volcano-related pollutants from Etna volcano, and succeeded in predicting their emission, transport and settling in the Mediterranean area. In the case of Etna volcano, different efforts have been made in the search of new techniques to model volcanic plumes using different tools, such as Embedded Computer Vision (ECV) camera, LiDAR measurements, and more precise estimations of ash concentrations with the intention of complementing VAAC’s efforts (Scollo et al. 2015). Continuous efforts are made to monitor volcanic plumes from Etna, and any new perspective that aids in tracking their trajectories greatly enhances these surveillance endeavours. This monitoring task is ongoing and essential (Toulouse VAAC 2021).

Governmental agencies such as the INGEMMET-OVI (Observatorio Volcanológico del Instituto Geológico, Minero y Metalúrgico; Peru INGEMMET, 2020) provide daily ash dispersal forecasts for some frequently active Peruvian volcanoes (e.g., Sabancaya) by using the Ash3d tool (Mastin et al. 2021). In this case, the different modelled eruptive scenarios are shown on their webpage.

When volcanoes are situated in areas with high population density and significant infrastructure presence (e.g., airports, power plants, hospitals) (Malawani et al. 2021; Alberico et al. 2011)., it is imperative to accurately forecast ash dispersion and deposition, such is the case for Popocatepetl Volcano near Mexico Megacity (e.g. De La Cruz-Reyna and Tilling 2008) around it large towns and villages are within 20 km from the crater with > 165,000 inhabitants; (INEGI, 2011). Ash dispersal is usually reported on the Washington VAAC webpage when concentration allows ash detection and tracking. Popocatépetl volcano, located ~ 60 km southeast from Mexico Megacity, has shown persistent eruptive activity since 1994 (e.g., Delgado-Granados et al. 2001; Martin-Del Pozzo et al. 2008). Popocatépetl’s current activity includes strong degassing (e.g., Matiella Novak et al. 2008; Grutter et al. 2008), ash emissions and lava eruptions, raising concern on the local authorities about the possible effects of volcanic ash emissions, particularly on the impact of ash fall on the population’s health. Another matter of concern is the public infrastructure; for example, the International Airport of Puebla (AIPU) 30 km from the crater of the volcano, Mexico City Airport (AICM) is 65 km away, Cuernavaca City Airport (ACV) is 70 km, and the new International Airport Felipe Ángeles (AIFA) within 90 km in distance and important highways that connect the north of the country with the south; Fig. 1). The national government still does not have a system for decision-making and for implement precautionary measures against ash cloud exposure with sufficient lead time for preventive action.

Fig. 1

Location map of Popocatépetl volcano, Mexico City (CDMX), Puebla City (CPUE), and Cuernavaca City (CDCV) (shaded in grey). Location of the most important airports surrounding Popocatépetl volcano are shown: The Airport of Cuernavaca City (ACV, 70 km away), the International Airport of Puebla (AIPU, 30 km away), the international Airport of Mexico City (AICM, 65 km away), and the International Airport Felipe Ángeles (AIFA, 90 km away)

Under a typical wind velocity of ∼10 m/s, volcanic ash may reach the closest villages in less than 20 min and the Puebla international airport in less than 1 h. Therefore, an important question to answer is, in which direction will the ash cloud move away from the volcano during an ash emission event? The area potentially affected by the volcano is very large with many infrastructure elements within 70 km from the crater (Macias et al. 1995; Martin Del Pozzo et al. 2017), so it is very difficult to apply preventive measures to all sectors at once due to the high cost of countermeasures and possible loss of credibility. In consequence, the National Centre for Disaster Prevention (CENAPRED), the official institution in charge of volcanic monitoring in Mexico, asked the Universidad Nacional Autónoma de México (UNAM) to develop a system aimed at forecasting ash dispersal under credible eruptive scenarios and realistic wind patterns.

This contribution describes the system developed and deployed at the facilities of CENAPRED.

Methods

To perform the ash dispersion forecast, a 36-h weather forecast has to be generated first. Subsequently, 108 scenarios are modelled corresponding to column heights of 3, 5 and 10 km above the crater starting every hour for the following 36 h. The selected heights are derived from observations of recent volcanic activity.(e.g., Delgado-Granados et al. 2001; Martin-Del Pozzo et al. 2008) and the expertise of civil protection personnel. Among the observed eruption phases, 3 km-high eruption columns are the most frequently recorded, while 10 km-high columns represent some of the highest observed during the current eruption phase. The 5 km-high columns fall between the other altitudes observed, being more likely to occur than the 10 km-high columns but less frequent than the 3 km-high columns. The duration of the plume emission for all these columns lasts approximately one hour. Each of these 108 scenarios is run for 8 h after the start time. In this way, a wide spectrum of eruption scenarios is considered.

Ash dispersion forecast depends on both wind conditions and eruption source parameters. In the following sections we describe the meteorological model, the ash dispersion and deposition model, the eruption source parameters, and the ash forecast system. Additional technical details are described in García-Reynoso et al. (2019).

Meteorological model

An important input for our forecasting model is the weather conditions, which have to be representative of the region and time of the year to simulate. We address this here by using the Weather Research and Forecasting (WRF) model system (Skamarock et al. 2008), and the Advanced Research WRF (ARW) is the dynamic core used from the two cores available in the WRF model (AWR and NMM) (Skamarock et al. 2019). The ash dispersion forecast system uses the output of the local weather forecast (GIOA Pronóstico Meteorológico 2017), having the advantages of reducing calculation time and working on the same computer system, the disadvantage being that both model resolution and geographic extension are fixed. The domains considered for the weather modelling are presented in Table 1. In this table, the larger domain is used to downscale the weather conditions and calculate the boundary conditions for the smaller domain. Here we consider the domain with the highest spatial resolution to obtain the best description of the weather conditions; this area includes the Mexico City International Airport (AICM) and Mexico City, both of which are located within an area of 70 km radius around the volcano (Fig. 1). The domain’s centre for as dispersion is located at 19.252°N and 98.586°W. The initial and boundary global forecast data come from the Global Forecast System (GFS) (2015) provided by the National Center for Environmental Prediction (NCEP). The weather forecast is set up to run for 72 h starting at 00Z time.

Table 1 WRF model domain definition. Cell and mesh dimensions

The WRF parameterizations used here are: microphysics from Morrison et al. (2015), an update from the Thompson et al. (2008) scheme Long-wave Radiation Rapid Radiative Transfer Model RRTM (Mlawer et al. 1997); short wave radiation Goddard scheme (Chou and Suarez 1994); surface layer ETA similarity (Janić 2001); Noah land surface model (Ek et al. 2003); planetary boundary layer, the Mello-Yamada-Janjic scheme (Mellor and Yamada 1982); and cumulus parameterization, the Kain-Fritsch scheme (Kain 2004).

Ash dispersion and deposition model

Here we used the FALL3D model (v.7.1) to simulate ash transport and deposition. This model utilizes a Finite Differences (FD) explicit scheme to solve the advection–diffusion-sedimentation (ADS) equations, and it is a three-dimensional Eulerian model designed for the transport and deposition of volcanic ashes. The model employs terrain-following coordinates to account for the effects of topography on the movement and settling of volcanic ash particles. A comprehensive description of FALL3D model is given in Costa et al. (2006; 2013) and Folch et al. (2009; 2016).

FALL3D has been widely utilized in the mapping of volcanic hazards at Popocatépetl volcano, with previous publications documenting its applications (Martin Del Pozzo et al. 2017).

FALL3D requires input data such Eruption Source Parameters (ESP), terrain elevation, crater location and Total Grain Size Distribution (TGSD); specific details for the Popocatepetl volcano can be found in García-Reynoso et al. (2019). Model outputs are total ground load, ash ground thickness, aerosol optical depth (AOD), particulate matter air concentrations (PM10, PM_2.5), and other related products.

The domain used for the FALL3D modelling has 100 by 100 squared cells, each cell is ~ 4 km in size and contains 21 vertical levels. In terms of model parameters, the horizontal turbulence is set to a constant value; the vertical turbulence varies with elevation (see Ulke 2000); the particle setting velocity model is based on the Ganser model (Ganser 1993), and the terrain elevation is obtained from the WRF output file. The parameters pertaining to the ash plume, the volcano, and the FALL3D parameterizations used in this study are described in and García Reynoso et al. (2018) and provided in more comprehensive detail in the Complementary material.

Eruption source parameters

Eruption characteristics are a key input for ash dispersion and deposition forecast. The problem in our case is the impossibility to forecast the starting time, duration and intensity of any future eruptive pulse. In order to resolve that, an analysis of all previous eruptions (e.g., Delgado Granados et al. 2008; Cross et al. 2012; Quezada-Reyes et al. 2013) was made here in order to identify the characteristics of the possible events that may occur. As explained above, previous events can be described as one hour of emission under three column heights: 3, 5 and 10 km above the crater. Modelled eruptions are set to be sourced from the crater’s lowest rim (5,160 m above mean sea level; total altitude of the volcano 5,452 m above mean sea level), which is located at 19.0231°N and 98.6246°W.

A key parameter in ash dispersion modelling is the Total Grain Size Distribution (TGSD), which is determined through field sampling. Linares López et al. conducted multiple studies (1998, 1999, 2004), while Linares López (2001) focused on the volcanic ashes emitted from Popocatépetl during the ongoing eruptive stage, spanning from 1994 to 2000 (51 events). The granulometric data obtained from these studies, ranging in size from -1ɸ to 4ɸ. The TGSD employed in the modelling is a Gaussian distribution derived from the granulometry data reported in the aforementioned research.

The eruptive column is set to a mushroom-like shape (i.e. Suzuki type; Suzuki 1983), and the mass flow rate (MFR) is estimated from column height and wind field using the Woodhouse approach (Woodhouse et al. 2013). All ash emission considers an eruption of one-hour duration.

The ash aggregation effects are based on the Cornell model (Cornell et al. 1983), and we do not consider a gravity current effect in our modelling.

To account for all potential scenarios within a day, the system was configured to run for 36 events, encompassing 24 h of the current day and an additional 12 h from the next day. The first event commences at 12Z on the current day, allowing 12 h for the weather model to achieve stability. Each eruption event has a duration of one hour, while the dispersion simulation spans 8 h (as ash emissions remain within the domain for up to 8 h). As part of this setup, three different eruptive columns are forecasted for each hour. Consequently, there are a total of 108 model runs per day.

Forecast system description

A set of Cshell scripts were developed to perform the sequence of steps to run both WRF and Fall3D. A flowchart is presented in Fig. 2, illustrating the three main activities of the script: 1) Weather forecast, which downloads the meteorological data, pre-processes the inputs, and runs the weather model; 2) Ash dispersion modelling, the ash dispersion modelling process involves a series of pre-processing steps carried out by the Fall3D system, which includes weather data processing (SetDBS), ash size distribution (SetGrn), and ash emissions estimation (SetSrc). These pre-processing steps generate the necessary inputs for ash dispersion modelling. The Fall3D system then proceeds to execute the dispersion process.; 3) Visualization step, involves processing the outcomes of ash dispersion and deposition to create images. These images are then uploaded to the web page, making the results easily accessible and viewable to users. Figure 3 describes an ash forecast procedure that utilizes weather forecasts generated by the IOA group and involves the use of the Fall3D model to generate 108 simulations (3 × 36) per day.

Fig. 2

Activities and auxiliary programs used to perform the ash dispersion forecast

Fig. 3

Ash forecast procedure using weather forecast produced by IOA group and generating 108 (3 × 36) simulations per day by using the Fall3D

Regarding the visualization products of the ash forecast, it is important to note that the approach we present here differs from the creation of a typical hazards map. Hazards maps are typically developed for civil protection planning to mitigate the impact of eruptions of various magnitudes. In contrast, our focus lies in a distinct scale of hazards mapping, offering day-by-day forecasts of hypothetical eruptions for real-time decision-making. Additionally, this system serves as a tool for civil aviation decision-making in near-real time, which traditional hazards maps do not provide. Consequently, the colour code employed in this system for forecasting ash clouds varies from the one used in "traditional" volcanic hazard maps.

Results

The system has been operational since December 2014. To compare model results with observations, seven different eruptions were selected. These eruptions occurred between 2015 and 2022, and the archived ash advisories from the Washington-VAAC website (Gallina and Street 2004) were utilized. For each year, an eruption event from Popocatepetl was randomly chosen, and its corresponding KML file was obtained (except for the year 2020, where no KML files were available). Subsequently, the output file from the forecast was identified using the date, and the forecast time corresponded to a forecast that started 3 h before the advisory time. The eruption events that were used for comparison are presented in Table 2; those eruptions had an average column height of 6,100 m asml. The comparison was made using similar Flight Level concentrations between the VAAC advisory report and our forecast. Figures 4, 5, 6, 7, 8, 9 and 10 show the ash advisory (polygon) from VAAC and the ash concentration from the model (round shapes). By analysing these figures, it becomes evident that, on the whole, the model effectively captures the overall dispersion direction of the ash cloud. Furthermore, it accurately reproduces the extent of the ash's influence. Discrepancies between the model and observation can be attributed to the distinct resolutions utilized in each approach.

Table 2 Eruptions dates used for comparison for year 2015 to 2022. Source: CENAPRED

Fig. 4

Eruption 14:51 UTC, 14^th June, 2015. Observed ash, red polygon according to the Washington-VAAC. Forecasted ash concentration, coloured area in mg/m³(see legend). VAAC image taken from the VAAC-Washington webpage

Fig. 5

Eruption 18:46 UTC, 20^th May, 2016. Observed ash, red polygon according to the Washington-VAAC. Forecasted ash concentration, coloured area in mg/m³ (see legend). VAAC image taken from the VAAC-Washington webpage

Fig. 6

Eruption 14:20 UTC, 31 July, 2017. Observed ash, red polygon according to the Washington-VAAC. Forecasted ash concentration, coloured area in mg/m³ (see legend). VAAC image taken from the VAAC-Washington webpage

Fig. 7

Eruption 07:16 UTC, 1^st March, 2018. Observed ash, red polygon according to the Washington-VAAC. Forecasted ash concentration, coloured area in mg/m³ (see legend). VAAC image taken from the VAAC-Washington webpage

Fig. 8

Eruption 19:49 UTC 14 August, 2019. Observed ash, red polygon according to the Washington-VAAC. Forecasted ash concentration, coloured area in mg/m³ (see legend). VAAC image taken from the VAAC-Washington webpage

Fig. 9

Eruption 06:23 UTC 24^th March, 2021. Observed ash, red polygon according to the Washington-VAAC. Forecasted ash concentration, coloured area in mg/m³ (see legend). VAAC image taken from the VAAC-Washington webpage

Fig. 10

Eruption 06:22 UTC 21^st October, 2022. Observed ash, red polygon according to the Washington-VAAC. Forecasted ash concentration, coloured area in mg/m³ (see legend). VAAC image taken from the VAAC-Washington webpage

Operational system

A local webpage was built (GIOA Dispersión de Cenizas 2015) to display ash dispersion-deposition forecast results (Fig. 11). It uses the following technologies: CMS Joomla 3.4.1, Apache 2.4.7, PHP 5.5.9–1 and MySQL 5.5.46. The web page contains four tabs: 1) Participants, 2) model information, 3) user manual and 4) previous site. It includes the following products: flight levels from 100 to 400 (hundreds of feet above sea-level), surface concentration (g/m³), ash thickness (mm), column mass (g/m²), load deposition (kg/m²), Aerosol Optical Depth (AOD), a panel with all the flight levels, and to visualize the vertical distribution of ash, we utilize longitude-latitude and S–N vertical planes. These planes present the ash airborne concentration using the crater as the central reference point (PM₁₀ and total). The products were selected to provide information useful for decision makers. For visualization purposes, the flight level 400 is used based on the maximum height of commercial flights, even though ash transport can reach the flight level 500.

Fig. 11

A sample view of the project’s web page. http://grupo-ioa.atmosfera.unam.mx/pronosticos/index.php/dispersion-de-cenizas Image taken from GIOA webpage

There are different ways to represent the results from modelling, including colour schemes and additional hazard information on a map (Brewer 1994; Brewer et al. 2013; Thompson et al. 2015). Throughout the project's development, there was significant and ongoing interaction between potential users and developers. This valuable feedback influenced the selection of the deterministic scenario map for each product display. Additionally, a colour scheme of blue-green–red-brown was adopted to represent ash concentrations, and a logarithmic scale was utilized to display certain products, such as surface concentrations.

Discussion

In the field of ash forecasting, several studies have been conducted to address different aspects of the process. Heffter and Stunder (1993), Collini et al. (2013), Arnold et al. (2014) and Hasegawa et al. (2015) used the characteristics of an ongoing eruption to forecast ash dispersion. Connor et al. (2001), Jenkins et al. (2015) and Uesawa et al. (2022) used historical records to estimate the expected distribution of ash deposition. On the other hand, Parra and Folch (2015) focused on forecasting a single eruption case. Additionally, specific monitoring and forecasting systems based on monitoring data have been developed for Etna, Italy (Scollo et al. 2009) and Sangay volcano, Ecuador (Bernard et al. 2022); these studies also provide online forecasts and hazard maps so they can be consulted by civil protection authorities. However, in these systems the information has a certain delay from the eruption event to the ash dispersion forecast, so an operator and a large-capacity computer system are required to obtain results in short time. In contrast, our methodology has a faster response time because whenever an eruption occurs, there are already results for a set of several possible eruptions available, thus reducing the time to issue the ash dispersal and deposition forecast. The method developed here has the advantage of providing more information in real time to facilitate the decision-making process.

Our forecast system was developed for operational purposes, to provide information for technical monitoring personnel and, because of that, future work is required to develop a set of products focused on likelihood and probabilities of ash dispersion to the general public (Connor et al. 2001; Doyle et al. 2014). A fundamental objective of this system is to ensure the delivery of resulting maps directly to decision makers through their cellular phones or tablets, enabling them to access the results effortlessly at any time and from any location. This accessibility is crucial for triggering the protocols they have already designed in response to the forecasted ash dispersion scenarios.

To assess the dispersion performance, we conducted evaluations over a period spanning from 2016 to 2022. During this evaluation period, several eruptions were examined, and a comparison was made by contrasting our model against the Washington VAAC reports. The findings indicated a noteworthy similarity in the ash trajectory and dispersion patterns between our model and the VAAC reports.

In addition, the forecast has fixed ash density and size distribution having influence in the sedimentation velocity, if the ash density is higher in the forecast than in the observed eruption, the concentration shown by the model will occur at lower flight level than the real one. On the other hand, we always have to have in mind that the column height of the forecast could be different from the real one, and that affects the plume direction due to a wind shear.

The model developed here proved its utility in several eruptive events in 2016. For instance, during the night of July 31^st, Popocatépetl erupted explosively at 23:24 Local Time, producing an ash cloud that reached Mexico City until the early morning of August 1^st. Several sectors of southern Mexico City were partially covered by a thin film of fine ash. In anticipation, at 22:46, July 31^st (i.e., ~ 45 min before the eruption), the forecast produced by our system was released to both authorities and volcanologists of the Advisory Scientific Committee (part of the National Civil Protection System, known as SINAPROC).

This event gained significant attention among the people of Mexico City, as reported by the media (Gómez-Flores et al. 2016; Navarro 2016). Although there were numerous eruptions during this period, only a few of them actually reached the City, making this particular case noteworthy for its impact on the urban area. The general wind pattern indicates that half of the year the winds blow towards the east, and the other half towards the west, i.e., to Mexico City. The westward wind pattern coincides with the rainy season, so it is more difficult to directly observe the ash cloud moving towards the city, highlighting another importance of our forecasting system. Volcanic explosions can be detected through seismograms, but if the surveillance system is blind to see the ash cloud, our near-real-time forecast helps prevent the consequences of the ash fall by giving timely advice and online information.

The system developed here has demonstrated to be an important tool for decision-making. It has helped preventing the impact of the ashes by providing both civil protection and aviation authorities timely information. Civil protection authorities alert the local officials at every municipality to be prepared for volcanic ash accumulated on roofs, and advice people about ash mitigation measures. The recommendations made by Doyle et al. (2014) about ways to communicate the results to the public, albeit useful, are unfortunately out of the scope of this project because all communication protocols are already set by the Ministry of the Interior, and must be followed by the Civil Defence and CENAPRED at the Ministry of Security and Citizen Protection (SSPC). It is important to mention that an important agreement among authorities and scientists is that all communications about scientific matters should be made by the authorities, not by scientists, this in order to have a single communication channel with the general public. This agreement was established in 1994, at the beginning of Popocatépetl’s current eruptive stage.

Conclusions

Here we presented an automated ash forecast dispersion system for Popocatépetl volcano, Mexico. Popocatépetl volcano is frequently erupting since 1994, and has had several effusive and explosive events so far. Based on the previous knowledge about the volcano, our system models 108 possible eruptive scenarios on a daily basis. The set of scenarios consist of a combination of three column heights, starting every hour, for the next 36 h, and with every eruption activity lasting one hour. In order to evaluate the dispersion performance, a 2016–2022 evaluation period was used. where some eruptions were compared contrasting our model against the Washington VAAC reports. The results showed that both model and VAAC reports have similar ash trajectory and dispersion. Our system has been successfully used by civil protection authorities, supporting decision-making.

The automated method presented generates a range of high-resolution modelling scenarios, encompassing a common set of potential eruptions. As a result, it eliminates the need for an operator with specialized modelling skills or extensive computing power to obtain rapid and dependable forecasts. Additionally, the system holds the advantage of easy adaptability to other similar cases.

Our system provides decision makers in at-risk areas with additional information, so civil protection alerts can be issued timely. Additional work is still required to generate products aimed at the general public.

This system can be adapted to any other active and hazardous volcano. For instance, the same system could be adapted to forecast the distribution of ash clouds from Colima volcano in western Mexico. This is a forthcoming step to follow.