Volcanic hazard assessment for tephra fallout in Martinique

Abstract

Mount Pelée (Martinique) is one of the most active volcanoes in the Lesser Antilles arc with more than 34 magmatic events in the last 24,000 years, including the deadliest eruption of the 20^th century. The current volcanic hazard map used in the civil security plan puts the emphasis on the volcanic hazard close to the volcano. This map is however based on an incomplete eruptive history and does not take into account the variability of the expected source conditions (mass eruption rate, total erupted mass, and grain-size distribution) or the wind effect on ash dispersal. We propose here to refine the volcanic hazard map for tephra fallout by using the 2-D model of ash dispersal HAZMAP. We first simulate the maximum expected eruptive scenario at Mount Pelée (i.e., the P3 eruption) using a seasonal wind profile. Building upon the good agreement with field data, we compute probability maps based on this maximum expected scenario, which show that tephra fallout hazard could threaten not only areas close to the volcano but also the southern part of Martinique. We then use a comprehensive approach based on 16 eruptive scenarios that include new field constraints obtained in the recent years on the past Plinian eruptions of Mount Pelée volcano. Each eruptive scenario considers different values of total erupted mass and mass eruption rate, and is characterized by a given probability of occurrence estimated from the refined eruptive history of the volcano. The 1979-2019 meteorological ERA-5 database is used to further take into account the daily variability of winds. These new probability maps show that the area of probable total destruction is wider when considering the 16 scenarios compared to the maximum expected scenario. The southern part of Martinique, although less threatened than when considering the maximum expected scenario, would still be impacted both by tephra fallout and by its high dependence on the water and electrical network carried from the northern part of the island. Finally, we show that key infrastructures in Martinique (such as the international airport) have a non-negligible probability of being impacted by a future Plinian eruption of the Mount Pelée. These results provide strong arguments for and will support significant and timely reconceiving of the emergency procedures as the local authorities have now placed Mount Pelée volcano on alert level yellow (vigilance) based on increased seismicity and tremor-type signals.

Introduction

Mount Pelée (Martinique) is one of the most active volcanoes in the Lesser Antilles arc with more than 34 magmatic events in the last 24,000 years (Smith and Roobol 1990; Westercamp and Traineau 1983; Boudon et al. 2005; Michaud-Dubuy 2019). Mount Pelée is mainly known for the 1902 −1905 Pelean eruption during which two pyroclastic flows were responsible for the complete destruction of Saint-Pierre and Morne-Rouge as well as for the deaths of nearly 30,000 people, making this eruption the deadliest of the 20^th century. The recent activity of Mount Pelée volcano was however not only characterized by Pelean-style eruptions but also by more powerful Plinian-style events, which impacted areas that are now densely populated (Michaud-Dubuy et al. 2019). Today, about 400,000 people live in Martinique and are thus threatened by volcanic hazards ranging from pyroclastic flows to distal tephra fallout, which makes crucial the improvement of volcanic hazard assessment (and consequently risk assessment) especially in a context where Mount Pelée exhibits signs of increased activity since 2020. Note that as in Bonadonna and Costa (2013), we use in this article the word “tephra” as a collective term for all particles ejected from volcanoes, irrespective of size, shape, and composition, whereas “tephra fall” refers to the process of particle fallout.

Stratigraphic studies show that Mount Pelée volcano is characterized by a heterogeneous succession of two eruptive styles in the last 25 ka: Pelean (dome-forming eruptions), and Plinian (open-vent eruptions producing a sustained eruptive column) (Roobol and Smith 1976; Westercamp and Traineau 1983). The Pelean events include 21 eruptions that generated pyroclastic flows mostly confined into the deep valleys around Mount Pelée and occasionally filled them entirely (e.g., the “rivière Blanche” in 1902, Lacroix (1904)). In between these major eruptions, Mount Pelée produced more powerful Plinian eruptions (Westercamp and Traineau 1983; Traineau et al. 1989; Bardintzeff et al. 1989) and underwent three major flank collapses (Le Friant et al. 2003; Solaro et al. 2020). The Plinian eruptive history of Mount Pelée was extensively revisited in the last decade (Carazzo et al. 2012; Carazzo et al. 2019; Carazzo et al. 2020; Michaud-Dubuy et al. 2019) in order to complete the recent eruptive history of the volcano (Michaud-Dubuy 2019). These studies, based on more than 200 outcrops and robust age determinations highlight that Plinian eruptions are more frequent at Mount Pelée than previously thought, and provide important Eruptive Source Parameters (ESPs, total volume, mass eruption rate, total grain-size distribution, column height, exit velocity) retrieved from field observations and measurements combined with physical models of volcanic columns. Another major outcome is that the associated hazards, including tephra fallout whose dispersal strongly depends on wind speed and direction, and pyroclastic flows due to column collapse, are largely underestimated in the current evacuation plans (Michaud-Dubuy 2019). Although the most probable eruptive scenario for Mount Pelée volcano in the future may be a phreatic eruption (Boudon et al. 2005), the possibility of a more powerful Plinian eruption occurring is non-negligible and could threaten the entire Caribbean region.

The characteristics of the wind have a primary control on the dispersion of volcanic tephra and thus on the associated hazards. Due to its central location in the Lesser Antilles arc, the island of Martinique (14 ^∘40"N, 61 ^∘00"W) is dominated by an oceanic tropical climate that can be divided into two main seasons (Gouirand et al. 2020): the summer rainy season extending from April/May to October/November (during which the cyclonic hazard is higher), and the winter dry season (also named Lent season) extending from November to April. The tropical Caribbean climate is characterized by a strong influence of the northern hemisphere easterly trade winds in the low to mid-troposphere (from the surface to ≈ 5 −7 km-high), and of the counter-trade winds in the upper troposphere (between 7 and 17 km-high). The high variability of these winds from the summer to the winter season, in both direction and speed, has often had a strong impact on tephra dispersal during an eruptive event (Michaud-Dubuy et al. 2019) as any sufficiently high volcanic column is affected by winds (Komorowski et al. 2008). Moreover, our recent study on the Bellefontaine Plinian eruption of Mount Pelée shows that specific wind conditions can lead to the dispersal of tephra to the south of the island, a direction that is not expected from the season-averaged or the monthly-averaged wind profiles. It is therefore crucial to account for the daily variability of winds to produce robust hazard maps in Martinique (Michaud-Dubuy et al. 2019).

During the last few decades, various approaches have been proposed to assess volcanic hazards and to improve the management of volcanic crises. A traditional approach to anticipate the impact of a future eruption relies on geological records to produce hazard maps (Baker 1985; Houghton et al. 1987; Stieltjes and Mirgon 1998; Newhall and Hoblitt 2002; Orsi et al. 2004). This method is often the only one possible for volcanic environments where only limited information is available on the past eruptive history of the volcano, but it can hide crucial information about the weakest (and more frequent) events, as their deposits are rapidly eroded or buried beneath those of more voluminous eruptions. Since the 1990’s, Geographic Information System (GIS) tools and numerical simulations can complement the geological approach by providing quantitative information. Numerical modeling, combined with field studies, can be used to simulate an historical eruption (Suzuki and Koyaguchi 2013; Gueugneau et al. 2020; Vicari et al. 2007), and to further explore a wide range of possible scenarios in order to identify potential affected areas (Bonadonna et al. 2002; Macías et al. 2008; Macedonio et al. 2016). Tephra fallout hazard assessment commonly relies on dispersal models fed by geological information retrieved in the field. Previous studies focused on either one eruptive scenario (usually the largest one or the more likely to happen within a given time window) along with a large set a wind profiles (Komorowski et al. 2008; Folch and Sulpizio 2010; Bonasia et al. 2011), or considered instead one or several eruptive scenarios together with a single wind profile that is commonly averaged over a season or estimated to be the most probable one (Barberi et al. 1992; Tsuji et al. 2017). The most widely used method consists of considering several eruptive scenarios as input parameters along with a wide set of wind profiles (Cioni et al. 2003; Bonadonna et al. 2005; Macedonio et al. 2008; Costa et al. 2009; Bonasia et al. 2012; Bonasia et al. 2014; Becerril et al. 2014; Macedonio et al. 2016; Biass et al. 2017). The outcomes of these studies generally take the form of hazard maps showing probabilities of reaching a tephra loading (in kgm⁻²) greater than a given threshold.

In this study, we aim at revising the current hazard map for tephra fallout at Mount Pelée volcano by accounting for its eruptive history and the specific wind conditions in the Lesser Antilles. For this, we use the 2-D dispersal model HAZMAP (version 2.4.2), together with several eruptive scenarios based on the stratigraphical record (and including the maximum expected scenario for Martinique) and the ERA5 dataset that accounts for the daily variability of winds. We chose to use the HAZMAP model for its low computational cost (a single simulation can be run in less than 5 seconds) and its ease of use. First, we describe the current hazard map for tephra fallout. We then simulate the maximum known eruptive scenario at Mount Pelée and generate probability maps based on this scenario. Then, we describe the other Plinian eruptions identified and reconstructed at Mount Pelée, which allow us to build a matrix of 16 eruptive scenarios. The latter are used, together with an extensive wind database, to produce robust probability hazard maps for tephra fallout. Finally, we discuss the consequences of tephra fallout on critical infrastructures such as the international airport of Fort-de-France, and the possible occurrence of other volcanic hazards in Martinique.

Current hazard map for tephra fallout

Westercamp (1983) first proposed a preliminary zonation of volcanic hazards associated with PDC and fallout deposits in Martinique, on the basis of field observations. The current volcanic hazard map used in the ORSEC plan (the national emergency response plan) was developed by Stieltjes and Mirgon (1998). This map is also included in the Volcanic Hazard Atlas of Lesser Antilles (Boudon et al. 2005) and in the Department Document on Major Risks registering all natural risks at the scale of the island (DDRM 2014). To produce it, Stieltjes and Mirgon (1998) mapped the hazard zoning of each volcanic phenomenon considered in Martinique (i.e., tephra fallout, pyroclastic flows, lava intrusions/flows, gas emissions, lahars, landslides, and tsunamis), by using “exposure” matrices. These matrices combine both the intensity (I) and the frequency (F) of each volcanic phenomenon over the entire area exposed to it. Five classes of intensity and frequency are proposed by Stieltjes and Mirgon (1998), from I0/F0 for the lowest one to I4/F4 for the highest one, based on the past eruptive history of the Mount Pelée volcano known at that time (i.e., the last 5,000 years; Westercamp and Traineau (1983)). Seven hazard zoning maps were created, one for each of the seven volcanic phenomenon considered, and combined into the final integrated volcanic hazard map.

As the goal of the present work is to re-assess the tephra fallout hazard in Martinique, we only describe in detail the current hazard map for tephra fallout produced by Stieltjes and Mirgon (1998) and presented in Fig. 1. Five classes of intensity (from I0, very low to non-existent, to I4, very high) and one class of frequency F are taken into account in this map; the color scale thus depends on the five levels of exposure to tephra fallout hazard defined by the product I × F (F being equal to 1). Figure 1 shows that the northern part of the island is the most exposed to tephra fallout, and that the exposure level decreases with distance from the Mount Pelée summit. The southern half of Martinique is considered to be safe, as the exposure level is null (white color) beyond the Lamentin plain (airport location).

Fig. 1

Current hazard map for tephra fallout in Martinique, based on data from BRGM and built by Stieltjes and Mirgon (1998). Colors correspond to the exposure level with red: very high; pink: high; orange: intermediate; yellow: low; and white: very low to null. All maps were generated using the open source QGIS software. Coordinates are in WGS 84 − UTM Zone 20N system

The methodology used by Stieltjes and Mirgon (1998) is subject to two limitations. First, it is not based on numerical modeling of volcanic column dynamics and tephra dispersion caused by winds, as apparent from the circular shape of the exposure level areas. Second, the methodology was built on the eruptive history determined by Westercamp and Traineau (1983) for the last 5,000 years, which was since revisited and completed by Carazzo et al. (2012); Carazzo et al. (2019); Carazzo et al. (2020); Michaud-Dubuy et al. (2019). Recent tephrostratigraphical studies indeed allowed to identify, interpret and reconstruct the dynamics of nine Plinian eruptions in the last 24,000 years. Four of them, the P1 (1,300 cal CE) (Carazzo et al. 2012), P2 (280 cal CE) (Carazzo et al. 2019), P3 (79 cal CE) (Carazzo et al. 2020), and P5 (4,534 cal BP) (Westercamp and Traineau 1983) eruptions, are interpreted to be sub-Plinian to Plinian eruptions with the formation of a 20 to 30 km-high stable plume that underwent total or partial collapse with the production of associated pyroclastic flows at some stages (Carazzo et al. 2012; Carazzo et al. 2019; Carazzo et al. 2020). The air fall deposits from these eruptions can be observed on the volcano flanks, separated from each other by soil deposits attesting that they are not from the same eruption (Fig. 2). They present a large variability in their dispersal direction: tephra from the P5 (4,534 cal BP) and P2 (280 cal CE) eruptions blanketed the eastern flanks of the volcano (Westercamp and Traineau 1983; Carazzo et al. 2019) while those of the P1 (1,300 cal CE) and P3 (79 cal CE) events were deposited mostly on the western flanks (Carazzo et al. 2012; Carazzo et al. 2020). The Bellefontaine (13,516 cal BP), Carbet (18,711 cal BP) and Etoile (21,450 cal BP) eruptions are much older, and characterized by a stable column that dispersed tephra on the southern flanks, which we interpret to be due to the occurrence of winds of more peculiar orientation (Michaud-Dubuy et al. 2019; Michaud-Dubuy 2019). The P6 (4,610 BP) and P4 (2,440 BP) Plinian eruptions produced small collapsing fountains associated with pyroclastic flows filling several western valleys (Westercamp and Traineau 1983). In this work where we seek to assess the tephra fallout hazard, we focus on Plinian eruptions that generated a stable column, and discard both Pelean eruptions and collapsing stages of Plinian events, for which the impact of tephra fallout is less important.

Fig. 2

Representative photographs of outcrops of pumice fallout deposits from recent Plinian eruptions in Martinique: the 1,300 cal CE P1 (Carazzo et al. 2012), 280 cal CE P2 (Carazzo et al. 2019), 79 cal CE P3 (Carazzo et al. 2020), and 4,534 cal BP P5 events. The outcrops are presented from the closest to the volcano (left, a, 4.7 km) to the furthest (right, c, 9 km). Outcrop b is located at a medium distance (5.5 km). All scale bars are 25 cm long. Boundaries between units are marked by white dashed lines

Maximum expected eruption: VEI5 P3 event

Reliable hazard assessment can only be achieved by first testing and calibrating the model used on a well-known eruption (Bonadonna and Costa 2013). In this work, we performed all simulations with HAZMAP-2.4.2 (Macedonio et al. 2005), a semi-analytical model solving the equations of dispersion, transport and sedimentation of tephra and now routinely used for volcanic hazard assessment (Macedonio et al. 2008; Macedonio et al. 2016; Komorowski et al. 2008; Capra et al. 2008; Costa et al. 2009; Bonasia et al. 2011; Bonasia et al. 2012; Becerril et al. 2014). In a previous study, we successfully reproduced the Bellefontaine eruption with this model (Michaud-Dubuy et al. 2019), validating the methodology for moderate Plinian eruptions in Martinique. In this section, we seek to test the model against the most voluminous and powerful eruption recorded in Martinique: the 79 cal CE P3 event (Carazzo et al. 2020). There is no evidence of a larger eruption during the last 20,000 years in the Lesser Antilles. Prior to that time, only a few VEI5 ignimbritic eruptions occurred in the central part of Dominica producing slightly more voluminous products (Howe et al. 2014; Boudon et al. 2017). We therefore consider the P3 eruption to be the most catastrophic scenario that may occur in Martinique. Two sets of input parameters are required to run the simulations: ESPs and wind fields. We first simulate the P3 eruption using HAZMAP in its deposit mode (and thus requiring a single wind profile). We then use the probabilistic mode to compute several probability maps based on this maximum expected scenario. This mode requires several wind profiles including wind velocity components (u,v) as a function of altitude. Wind fields can either be taken from radiosonde measurements or from reanalysis data, which is the method we chose for this study. Finally, HAZMAP does not take into account particle aggregation, a phenomenon that we did not observe in the field.

Simulation of the P3 eruption

The P3 event can be divided into seven major Plinian phases that produced a total of 1 km³ dense rock equivalent (DRE) of deposits, which makes this event a VEI5 eruption (Carazzo et al. 2020). The first phase, named P3A, produced a stable column that covered with tephra the western slopes. The eruption then evolved towards a more unstable regime of partial column collapse generating pyroclastic flows. The fall deposits (phases P3A-C-E-G) can be found on both the western and eastern sides of the volcano whereas the pyroclastic flow deposits (P3B-D-F) are confined in western valleys (Carazzo et al. 2020). We seek to reproduce the main fall deposit of this event, P3A, with the HAZMAP model in its deposit mode (i.e., determinist and only requiring a single wind profile). The model requires several ESPs to run the simulations. The deposit mass, the maximum column height and the total grain-size distribution (TGSD), the three key parameters for tephra dispersal, were retrieved from field data and published in a previous study (Carazzo et al. 2020). We summarize below these volcanological parameters, together with the other inputs required by HAZMAP and used in our simulation. More detail about the methods used to combine our field observations and physical models of a volcanic plume in order to estimate the ESPs and their uncertainties can be found in Carazzo et al. (2020).

The maximum column height estimated from the distribution of lithic fragments by Carazzo et al. (2020) is 30 km. The volume of the P3A unit inferred from the thinning behavior of the deposits together with best fit functions (Daggitt et al. 2014) is 0.1 km³ DRE (Carazzo et al. 2020). We take a total mass of 2 ×10¹¹ kg, assuming a deposit density of 1,070 kgm⁻³ based on previous estimates for Mount Pelée deposits (Traineau et al. 1989) and consistent with the densities used in the literature (between 900 and 1500 kgm⁻³, Spence et al. (2008); Bonasia et al. (2011); Bonasia et al. (2012)). The TGSD was reconstructed by Carazzo et al. (2020) using deposit samples from several outcrops located all around the volcano together with the method of Kaminski and Jaupart (1998). The TGSD was found to be bimodal with a primary fine mode at 2 ϕ (ϕ being a particle size notation with the particle diameter in millimeters d_ϕ = 2^−ϕ) and a secondary coarse mode at -2 ϕ, and the median diameter and sorting are -0.2 ϕ and 2.3, respectively (Table S1). We note that the TGSD reconstructed from samples located on land is depleted in very fine ash material (ϕ>4), suggesting that most of the finest material was dispersed in the sea. The other input parameters required by HAZMAP are maintained constant in all the simulations presented in this paper. The mass distribution of particle in the volcanic column is parameterized using two Suzuki parameters that we set at A = 4 and λ = 1, as they represent a ratio H_B/ H_T (where H_B is the neutral buoyancy height of the plume, and H_T its maximum height) similar to the one observed for buoyant plumes (Morton et al. 1956; Sparks 1986). This set of Suzuki parameters is therefore the most commonly used in the literature as it also correctly reproduce the observed tephra deposits (Macedonio et al. 2008; Costa et al. 2009; Bonasia et al. 2011). In all simulations we use a horizontal atmospheric diffusion coefficient of 3000 m²s⁻¹, as this is the best-fit parameter found for simulating the Bellefontaine eruption in Michaud-Dubuy et al. (2019).

We first simulate the P3 eruption with HAZMAP in its deposit mode using the volcanological parameters estimated by Carazzo et al. (2020) and a typical wind profile in the Lesser Antilles. The main direction of dispersion suggests that P3 occurred during the rainy season. We thus use the mean rainy season wind profile in Michaud-Dubuy et al. (2019) characterized by north-westerlies with an average azimuth of 305 ^∘ in the upper troposphere (Fig. 3a). The results obtained are presented in Fig. 3b and c. Figure 3b reveals a good agreement between ground loads computed by HAZMAP and ground loads observed in the field by Carazzo et al. (2020). The simulated isopachs show a westward dispersal axis (Fig. 3c) that compares well with the one determined by Carazzo et al. (2020) for the P3A phase (their Fig. 6b). This good agreement confirms that HAZMAP can be confidently used to simulate tephra dispersion during large Plinian eruptions, and can thus be used to assess the volcanic hazard for tephra fallout when considering the maximum expected scenario at Mount Pelée volcano.

Fig. 3

Simulation of the P3A eruption (79 cal CE) with the HAZMAP model using the eruptive parameters given in the main text, and a a seasonal average wind profile (rainy season from Michaud-Dubuy et al. (2019), see main text). The results are presented with b the log-log plot of the observed versus calculated tephra loads (kgm⁻²) and c the isopach map obtained from the HAZMAP simulation with computed deposit thickness indicated in centimeters. Dashed lines in b indicate over- or under-estimations of 1/3 and 3 times the observed values, respectively

Meteorological dataset

In order to produce probability hazard maps for Mount Pelée volcano, we use wind azimuth and velocity profiles extracted from the European Centre for Medium-Range Weather Forecasts ERA-5 reanalysis for the years 1979-2019 (Hersbach et al. 2019). The initial content of ERA5 files consists of hourly global fields of zonal and meridional winds at a horizontal resolution of 0.25 ^∘×0.25 ^∘ (≈ 31 km) and vertically distributed on 37 pressure levels from 110 m (1000 hPa) to ≈ 48 km (1 hPa). These wind fields have been interpolated to match HAZMAP format by converting each of the 37 pressure levels into an altitude level using the altitude model (Figure S1). We select the wind components over Martinique at each time step and each pressure level in an area ranging from 14.7 ^∘N to 14.8 ^∘N and from 61.1 ^∘W to 61.2 ^∘W, and we calculate a mean daily wind profile for every day from January 1, 1979 to December 31, 2019. Our final dataset is thus composed of 14,975 vertical wind profiles (365 or 366 days times 41 years) shown in Fig. 4.

Fig. 4

Compass roses representing the 41-year wind database from the European Centre for Medium-Range Weather Forecasts ERA5 reanalysis (Hersbach et al. 2019) for a, c, e, the winter season and b, d, f, the rainy season. Horizontal wind vectors (speed and azimuth) are reported at different intervals of height: 0 to 7 km (a and b), 7 to 17 km (c and d), and 17 to 50 km (e and f). The wind speed is discretized into several layers, from 0 ms⁻¹ at the center of the rose to the maximum value of each dataset at the rose boundary

The lower troposphere (from the surface to ≈ 7 km in altitude) is characterized by E −NE trade winds during the winter season, with a maximum speed that is generally comprised between 0 and 15 ms⁻¹ but that can reach 27 ms⁻¹ in some rare cases (Fig. 4a). The upper troposphere (from 7 to 17 km) is characterized by strong westerlies (i.e., originating from the west) counter-trade winds with a speed reaching 40 ms⁻¹, thus of higher mean velocity compared to the trade winds (Fig. 4c). In the stratosphere (>17 km, Fig. 4e), the wind often comes from the east or the west, with a strong variability in speed (up to 80 ms⁻¹ in some rare cases). These stratospheric variations do not strongly influence tephra dispersal (and thus hazard assessment), as demonstrated in Michaud-Dubuy et al. (2019).

During the rainy season, the lower troposphere (up to ≈ 5 km) is characterized by constant E −SE trade winds, with speeds lower than for the winter season and generally ranging from 0 to 10 ms⁻¹ (Fig. 4b). In the upper trophosphere, winds mostly come from the west (even if it seems less pronounced than during the winter season) with a maximum wind speed often reaching 20 ms⁻¹ (Fig. 4d). Both wind speed and direction strongly vary in the stratosphere throughout the rainy season, with an azimuth of 90 ^∘ or 270 ^∘ and a maximum wind speed up to 80 ms⁻¹ (Fig. 4f). In the following parts, we use these wind profiles to assess the volcanic hazard for tephra fallout in Martinique.

Probability map with daily winds

We now use the model HAZMAP in its probability mode to compute ash loading probability maps for the maximum expected eruption scenario (i.e., the P3 event). We use the same volcanological parameters as for the P3 simulation, together with the complete 41-year wind data database (i.e., containing 14,975 daily wind profiles, Fig. 4). HAZMAP used in its probabilistic mode also requires static load threshold values to compute output hazard maps. We chose to use five thresholds based on the literature, each corresponding to a tephra thickness and a degree of damage on vegetation or infrastructures (Table 1). According to Wilson et al. (2017), a 1 mm-thick ash deposit (thus corresponding to a minimum loading of 1.07 kgm⁻²) leads to the airport closure, and to a required maintenance on all kinds of supply networks (e.g. electricity, water, wastewater, roads) in order to prevent further damage. At this stage, the visibility is already reduced. The second threshold, set at 10.7 kgm⁻² (i.e., 1 cm thickness), corresponds to extensive repair needed on supply networks as water for example is likely to be contaminated (Wilson et al. 2017). Other consequences are first damage to vegetation (Bonadonna et al. 2005), breathing difficulties, and car speed reduced by half. When the third threshold of 107 kgm⁻² (i.e., 10 cm thickness) is reached, replacement is required on supply networks and the airport is completely buried (Wilson et al. 2017). The main roads also become impassable for some vehicles. Beyond the fourth threshold of 214 kgm⁻² (i.e., 20 cm thickness), roads are impassable for all vehicles (Wilson et al. 2017), and roofs made of timber collapse (Komorowski et al. 2008; Baxter and Horwell 2015). The last threshold of 1,070 kgm⁻² (i.e., 1 m thickness) is the most critical as it corresponds to a complete destruction of all infrastructures and buildings (Wilson et al. 2017; Komorowski et al. 2008). The results for the P3 eruption are shown in Fig. 5.

Fig. 5

Ash loading probability maps for the maximum expected eruption scenario (P3 event) using the 41-year wind database, with a 5% probabilities of reaching the ash load thresholds of 10.7, 107, 214 and 1,070 kgm⁻², and b all probabilities to reach the 10.7 kgm⁻² threshold

Table 1 Tephra load thresholds considered for the hazard analysis, with their tephra thickness thresholds and corresponding damages on infrastructures, adapted from Bonadonna et al. (2005); Komorowski et al. (2008) and Wilson et al. (2017)

Figure 5a shows the 5% probability of reaching an ash loading corresponding to a thickness of 1 cm, 10 cm, 20 cm, and 1 m. The 5% probability of reaching a thickness of 1 mm is far south and is therefore not visible on this map. As a benchmark, the 5% probability of reaching a 1 cm-thick deposit corresponds to the ≈ 15 −20% probability of reaching a 1 mm-thick deposit. This map shows that the south of Martinique has a significant probability of being impacted by a future powerful eruption of Mount Pelée, a hazard that is not considered in the current hazard map (Fig. 1). The probability to reach a 1 cm thickness of ash is relatively low, but if that were to happen, some maintenance would still be required on various supply networks, especially on the water network as most of the drinkable water comes from the north of the island. The north of the island is the most threatened area because of the proximity with the volcano: there is indeed a non-negligible probability of reaching a thickness of 20 cm from Le Carbet (south of Saint-Pierre, Caribbean side) to Basse Pointe on the Atlantic coast (Fig. 5a). This thickness (corresponding to an ash mass load of 214 kgm⁻²) is critical as it corresponds to roof collapse of low- and medium-quality buildings. This risk is all the more important as a roof collapse and the consequent ash inflow may cause direct impact injury to the skull or body, suffocation and/or partial burial of the unhabitants (Baxter and Horwell 2015). Complete destruction may moreover be possible beyond Saint-Pierre, close to the volcano flanks, as there is a 5% probability of reaching 1 m thickness of ash.

Even a few millimeters of ash can cause major disruption, especially for the airport that is often the most effective way of evacuation in insular context. Figure 5b shows that there is a 15% probability of reaching 1 cm thickness of tephra in Fort-de-France and at the airport (and a 31% probability of reaching a tephra thickness of 1 mm), which means that a powerful eruption like the P3 event would most likely result in the airport closure. These conclusions are obtained for a VEI5 Plinian eruption, which occurred only once in the recent eruptive history of Mount Pelée (i.e., in the last 24,000 years).

The existence of a few explosive eruptions (∼63−30 ka Grand Bay, >22 ka Grande Savane, ∼80−51 ka Layou, ∼33−26 ka Roseau, and ∼30−24 ka Grand Fond, Howe et al. (2014); Boudon et al. (2017)) that produced a large volume of magma (2.5 −4 km³ DRE/eruption, Boudon et al. (2017)) in the last tens of thousands years in Dominica raises the question of the consequences of such a catastrophic scenario in Martinique. Figure S2 shows the 5% probability of reaching the ash load thresholds of 107, 214, and 1,070 kgm⁻² when considering a deposit mass of 10¹³ kg and a maximum column height of 39 km. With these extreme conditions, the model predicts at least 20 cm of deposits over the entire island of Martinique, and significant ash deposition in Dominica to the North and St. Lucia to the South. However, this scenario is very unlikely at Mount Pelée volcano due to the shallow crustal reservoir frequently drained by VEI3 −4 eruptions (Martel et al. 1998) preventing long magma accumulation timescales in reservoirs. The occurrence of very powerful explosive eruptions in Dominica may be related to specific structural and tectonic conditions creating extensive stresses that enhance the storage of large magma volumes (Jellinek and DePaolo 2003) at shallow (Howe et al. 2015) and/or deep crustal depths (Boudon et al. 2017). In any case, if such an eruption happened in the future in Dominica, Figure S2 shows that an inter-island cooperation would be of utmost importance for crisis management.

The next step is thus to assess the volcanic hazard for tephra fallout for other scenarios that are more likely to happen in the future. These scenarios are built from the geological record of VEI3-4 Plinian eruptions in Martinique, that we detail in the next section.

Geological record of VEI 3-4 Plinian eruptions in Martinique

Over the last 10 years, our team performed several extensive field studies in Martinique in order to reevaluate the eruptive history of the Mount Pelée and estimate key ESPs (Carazzo et al. 2012; Carazzo et al. 2019; Michaud-Dubuy et al. 2019; Michaud-Dubuy 2019; Carazzo et al. 2020). The resulting refined eruptive history is very rich, with the identification of 34 magmatic eruptions in the last 24 ka cal BP including 13 Plinian eruptions (Michaud-Dubuy 2019). The deposits of at least 6 Plinian eruptions were well-preserved enough to retrieve their eruptive parameters: the 1300 cal CE P1 (Carazzo et al. 2012), the 280 cal CE P2 (Carazzo et al. 2019), the 79 cal CE P3 (Carazzo et al. 2020), the 13,516 cal BP Bellefontaine (Michaud-Dubuy et al. 2019), the 18,711 cal BP Carbet (Michaud-Dubuy 2019), and the 21,450 cal BP Etoile eruptions (Michaud-Dubuy 2019). We summarize here the range of eruptive parameters estimated for these six Plinian events that we later use to define several eruptive scenarios for volcanic hazard assessment.

Eruptive parameters

All eruptive parameters and their uncertainties were retrieved from field data collected in Martinique (more than 200 outcrops studied on the flanks of Mount Pelée, see Michaud-Dubuy et al. (2019); Carazzo et al. (2020)) and published in Carazzo et al. (2012); Carazzo et al. (2019); Carazzo et al. (2020); Michaud-Dubuy et al. (2019); Michaud-Dubuy (2019). We do not use these estimates to simulate the impacts associated with a past eruption of Mount Pelée volcano, but to determine a typical range of values for each parameter that will be used, in turn, to define a scenario.

Total erupted mass

The volumes of the six studied Plinian eruptions are relatively close to each other comprising between 0.04 km³ DRE (for the Etoile eruption) and 1.02 km³ DRE (for the P3 eruption), which correspond to erupted masses ranging from ≈10¹¹ to ≈10¹² kg (Michaud-Dubuy 2019; Carazzo et al. 2020). All eruptions are thus ranked 4 on the VEI scale, with the exception of the P3 event, which is a VEI5 eruption (and the most powerful recorded at the Mount Pelée volcano). Bearing in mind that less powerful VEI3 eruptions most certainly occurred at Mount Pelée volcano, but did not leave any trace in the field due to severe erosion processes, we also consider lower masses (down to 10¹⁰ kg) in our hazard assessment.

Maximum column height and MER

The maximum column height was relatively similar for the 6 events with 19 −22 km for P1, 22 −26 km for P2, 19 −21 km for Bellefontaine, 19.6 km for Carbet, and 19 km for Etoile. The mass eruption rates (MER) associated with these column heights were also similar (≈10⁷kgs⁻¹). The P3 eruption again steps out from the usual pattern, with a maximum column height of 30 km and a MER >10⁸kgs⁻¹. As for the erupted mass, we also consider lower mass eruption rates (down to 10⁶kgs⁻¹) for the scenarios likely to happen in the future.

Total grain-size distribution

The total grain-size distributions of the 6 Plinian events were calculated using the method of Kaminski and Jaupart (1998), which accounts for the power-law size distribution of the rock fragments (i.e., where the numbers of particles with a radius larger than r is proportional to r^−D). The TGSD is then fully characterized by a power-law exponent D representing the fragmentation efficiency (details are given in Michaud-Dubuy et al. (2019)). It generally ranges from 2.9 (coarsest distribution) to 3.9 (finest distribution) for both fall and PDC deposits (Kaminski and Jaupart 1998). Grain-size analyses show that the eruptive products of the Bellefontaine are relatively coarser (power-law exponent D = 3.0 for the main unit) than those of the Carbet (D = 3.3), Etoile (D = 3.5), P3 (D = 3.3), P2 (D = 3.4 −3.5) and P1 eruptions (D = 3.2 −3.3). Because of their good preservation, we choose to only use the total grain-size distribution of the P1A (D = 3.2), P2A (D = 3.5), and P3A (D = 3.3) events (stable phases of the most recent eruptions) to calculate a mean TGSD used in all simulations (Table S1). This average TGSD is bimodal with a maximum peak at -3 ϕ and a secondary one at 1 ϕ. This average distribution allows to retain the fine grain-sizes preserved only for the P3 eruption (while mostly lost at sea for both P1 and P2 eruptions).

Definition of eruptive scenarios

Several methods can be used to construct volcanic hazard maps. For example, they can be based on the maximum expected eruptive scenario with a single or multiple wind profiles (Folch and Sulpizio 2010; Bonasia et al. 2011; Tsuji et al. 2017), two or three scenarios that are representative of the eruptive history of the volcano (Macedonio et al. 2008; Becerril et al. 2014), or a very large number of scenarios whose parameters are randomly sampled from a uniform or Gaussian probability distribution function (Bonadonna et al. 2005; Bonasia et al. 2014; Biass et al. 2017). In a previous section, we already assessed the volcanic hazard for tephra fallout when considering the maximum expected eruption at Mount Pelée volcano, based on the stratigraphic record (i.e., the VEI 5 P3 eruption). We now seek to build a new hazard map that accounts for the diversity of the possible eruptions at Mount Pelée volcano. This work is being carried out during the unrest of the volcano, which has been increased from alert level green (no alert) to level yellow (vigilance) in December 2020, based on the detection of an increased seismicity and tremor-type signals (OVSM-IPGP 2020). This change in the volcanic activity is giving us little time to characterize a full suite of eruption scenarios with appropriate distributions and parameterization, and we thus chose to consider a limited number of plausible and representative eruptive scenarios. By defining a larger set of scenarios, two major issues may arise in the resulting maps. The chosen set would consist of either very similar eruptive scenarios(and thus very similar tephra load maps), which may lead to a false degree of confidence in the final hazard map, or eruptive scenarios very different from each other resulting in a final hazard map with little practical value for the authorities (Jenkins et al. 2015). The final output map will be included in the revised version of the French emergency plan (ORSEC) in response to a volcanic disaster in Martinique, upon the request of the Prefecture of Martinique. Thus, the newly created volcanic hazard map must allow a realistic representation of the possible tephra mass loads on the island, especially as it will serve as a basis for the risk management and emergency evacuation plan updated by the Environment, Planning and Housing Agency (DEAL Martinique). To this aim, we define 16 eruptive scenarios using the method described below, and estimate their probability of occurrence based on the eruptive record of Mount Pelée.

Based on our refined eruptive history of the volcano, we conclude that the most probable future Plinian eruptive scenario in Martinique would be characterized by a MER comprised between 10⁶ and 10⁸kgs⁻¹, and a total erupted mass ranging from 10¹⁰ to 10¹² kg. We organize these MER and total mass ranges into a 4 ×4 matrix of 16 eruptive scenarios (Table 2). An eruptive scenario is defined in the matrix by a MER/Mass couple.

Table 2 Matrix of correlation used for the HAZMAP simulations, showing the relative probabilities of occurence of each eruptive scenario

We then estimate the probability of each eruptive scenario based on several assumptions:

• The probability of scenarios characterized by a total mass between 10¹¹ and 10¹² kg is twice the probability of cases with a total mass between 10¹⁰ and 10¹¹ kg, based on our stratigraphical record.

• Following the general observation that MER and total mass of deposits are positively correlated in Plinian eruptions (Carey and Sigurdsson 1989), we set a lower probability for scenarios with high MER/low total mass and low MER/high total mass, compared to scenarios characterized by a simultaneous increase or decrease in MER and total mass.

• Finally, we set a higher probability for the scenarios that are closer to the P1/P2/Bellefontaine eruptions (scenarios 10 and 11 in Table 2) as they represent the most frequent Plinian eruptive scenario at Mount Pelée volcano.

The calculated probabilities of occurrence for each of the 16 eruptive scenarios are reported in Table 2. Finally, as HAZMAP requires a maximum column height, we use our 1-D model PPM of volcanic column (Michaud-Dubuy et al. 2018) −considering a tropical atmosphere and an initial gas content of 3 wt %− to convert the MER values shown in Table 2 into maximum column heights. The gas content controls the column velocity at the vent and thus its stability, but it has no impact on the column height. The latter parameter varies between 13 and 24 km, a range of values consistent the maximum column heights retrieved from field data at Mount Pelée (except for the P3 eruption which is discussed separately, see section P3).

Figure 6 compares the probability distribution of our 16 eruptive scenarios to the one inferred from two worldwide explosive eruption databases: LaMEVE (Brown et al. 2014) and IVESPA (Aubry et al. 2021). We identified 135 events with information on both the MER and total erupted mass that can be used to sort them according to our 16 eruptive scenarios (Table S2). The resulting probability distribution is found to be in relatively good agreement with the one we built based on the eruptive record of Mount Pelée, which reinforces the confidence in the probabilities we set for each scenario.

Fig. 6

Probability distribution of the 16 eruptive scenarios based on the geological record at Mount Pelée volcano (in purple, see Table 2), and from two worldwide explosive eruption databases (LaMEVE and IVESPA, in grey; see Table S2)

Hazard maps for tephra fallout

Probability maps of selected scenarios

We performed 239,600 simulations (16 scenarios ×14,975 wind profiles) with HAZMAP in its probabilistic mode using our matrix of correlation (Table 2) and the daily wind profiles (Fig. 4). The Suzuki parameters (i.e., A = 4 and λ = 1) and horizontal atmospheric diffusion coefficient (i.e., 3000 m²s⁻¹) are the same than those used in the P3 simulation section. The results are presented using the five static load thresholds described in Table 1. The four maps in Fig. 7 show the results for the scenarios 4, 7, 10, and 13 characterized by a couple MER/Mass from the less powerful (scenario 4) to the most powerful (scenario 13). Only the 5% probability of reaching different static load thresholds is shown, and the missing thresholds are either largely exceeded on land, or not reached at all.

Fig. 7

5% probability maps for selected scenarios using the 41-year wind database, with a scenario 4, b scenario 7, c scenario 10, and d scenario 13 (Table 2). Different static load thresholds are shown for each scenario, note that missing thresholds are either largely exceeded on land or not reached at all

Figure 7 shows that the obtained hazard map strongly depends on the chosen scenario. From the minimum scenario considered (i.e., scenario 4, Fig. 7a) to the maximum one (i.e., scenario 13, Fig. 7d), the hazard level (and thus, the associated risks) are considerably different. Whereas only the northern part of Martinique has a 5% probability to reach a deposit thickness of 1 cm (i.e., an ash load threshold of 10.7 kgm⁻²) in the case of scenario 4 (Fig. 7a), this probability affect the entire island when considering the strongest scenario (Fig. 7d). This threshold (critical for the vegetation for instance) could affect half of the island for the scenario 7 (Fig. 7b), and could threaten the entire island down to Le Diamant in the case of scenario 10 (Fig. 7c). Note that the latter scenario (Fig. 7c) corresponds to the most probable one in the future, as the eruptive parameters considered resemble those of the P1, P2 and Bellefontaine eruptions. In this case, one can expect 5 cm of tephra at the airport, and more than 10 cm (i.e., 107 kgm⁻²) in the northern part of the island, which would cause potential contamination of water for the entire island. From the south of Saint-Pierre, the roads could be impassable for all vehicles because of the non-negligible probability to reach 20 cm of ash (i.e., 214 kgm⁻²) close to the volcano. When considering the strongest scenario (Fig. 7d), the entire northern part of the island could be completely destroyed (5% probability to reach 1 m of ash), and the international airport could be entirely buried under 10 cm of ash (i.e., 107 kgm⁻²). Note that even in the case of the weakest scenario (Fig. 7a), the Le Prêcheur area (on western flanks of the volcano) could be threatened by up to 10 cm of ash.

New hazard maps for Martinique

We now combine the simulation outputs for the 16 scenarios into a single map using the probabilities of occurrence of each scenario given in Table 2. The results are presented in Figs. 8 and 9.

Fig. 8

Ash load probability maps for Martinique using the 41-year wind database and considering the 16 eruptive scenarios, with probabilities to exceed a 1 cm of ash (10.7 kgm⁻²), b 10 cm of ash (107 kgm⁻²), c 20 cm of ash (214 kgm⁻²), and d 1 m of ash (1,070 kgm⁻²)

Fig. 9

Main hazard map for tephra fallout in Martinique using the 41-year wind database and considering the 16 eruptive scenarios, with a the 5% probabilities to reach a deposit thickness of 1 cm (10.7 kgm⁻²), 10 cm (107 kgm⁻²), 20 cm (214 kgm⁻²), and 1 m (1,070 kgm⁻²). Note that the 1 mm threshold (1.07 kgm⁻²) is largely exceeded on land and is therefore not visible. In b, the same map is represented using five levels of damage intensity with red: complete destruction, pink: very high, dark orange: high, pale orange: intermediate, and yellow: low. These levels depend on tephra load thresholds presented in Table 1

Figure 8 shows the probabilities of exceeding each static load threshold, except for the smallest one (i.e., 1.07 kgm⁻², corresponding to a 1 mm thickness of ash) as this threshold is largely exceeded in many locations of the island. As in the current hazard map for tephra fallout (Fig. 1), the northern part of Martinique remains the most threatened area in case of a future Plinian eruption at Mount Pelée volcano, with a 55% probability of reaching 1 cm of tephra (10.7 kgm⁻², Fig. 8a) at Saint-Pierre, a 25% chance of reaching 10 cm of tephra (107 kgm⁻², Fig. 8b), and even a 20% probability of exceeding 20 cm (214 kgm⁻², Fig. 8c) of ash. There is however less than a 5% probability of reaching the 1 m threshold (1,070 kgm⁻², Fig. 8d) at Saint-Pierre. These high probabilities however encompass mainly the Caribbean coast, as the winds in high altitudes (i.e., where most of tephra are injected) often come from the east. Basse Pointe, located on the Atlantic coast, has thus a lower probability to reach each threshold than Saint-Pierre, although it cannot be considered as a safe area from tephra hazard. Our stratigraphical record indeed shows that some Plinian eruptions (P2 and P5) mainly affected the eastern coast. Only the volcano flanks are likely to exceed the last threshold of 1,070 kgm⁻² and thus to be completely destroyed (Fig. 8d). This result is consistent with our field studies as we often measured more than 1 m of tephra in the northern part of Martinique (on every flanks of the volcano) (Carazzo et al. 2012; Carazzo et al. 2019; Michaud-Dubuy et al. 2019; Carazzo et al. 2020).

It is also important to note that the airport (and Fort-de-France, the most populated area in Martinique) exceeds a 10% probability of reaching the threshold corresponding to a tephra thickness of 1 cm (Fig. 8a), and thus has a even higher probability to be closed shortly after the beginning of a Plinian eruption. This event would have strong consequences on crisis management and evacuation procedures. The airport area has however less than 5% chance of being completely buried, and thus to be irremediably non-operational (Fig. 8b).

Finally, even if the southern part of the island seems to be safe from tephra fallout, it is still possible that tephra fallout occurs near Sainte-Anne depending on the dominant winds during the eruption. Our stratigraphical record indeed shows that at least one Plinian eruption occurred under northerly winds (Michaud-Dubuy et al. 2019) that affected the south of Martinique. In any case, since most of the drinkable water is carried from the northern part of the island, the southern part would still be largely affected, albeit indirectly.

Figure 9a shows the 5% probability of exceeding each threshold, and points out all the areas likely to be impacted by a future Plinian eruption of Mount Pelée. Note that the first threshold (1.07 kgm⁻²) is far exceeded and thus not visible on this map. The comparison between this map and the one obtained for the maximum expected eruption scenario (Fig. 5a) reveals two major differences. The area of probable total destruction (i.e., 1 m of tephra) is reduced when considering the maximum expected eruption scenario (∼77.5km²) compared to the 16 scenarios (∼190km²). On the contrary, the southern part of Martinique appears to be more threatened when considering the maximum expected scenario than when considering the 16 scenarios, as the 5% probability to reach a 1 cm tephra thickness encompasses Sainte-Anne in Fig. 5a but not in Fig. 9a. These differences are mainly due to the different TGSDs used in the simulations, as the P3 eruption is characterized by an overall finer TGSD than the one used for the simulations considering the 16 scenarios (Table S1). Figure 9b shows the expected damage intensity in Martinique using a color scale adapted from the intensity levels of Stieltjes and Mirgon (1998) to ease the comparison with the previous hazard map in Fig. 1. These results clearly demonstrate that the tephra fallout hazard needs to be considered in every part of Martinique.

Discussion

Tephra fallout on critical infrastructures

Disruption or damage to critical infrastructures can cause significant societal impacts and economic losses (Wilson et al. 2017). This is all the more critical in insular contexts such as Martinique, which is highly dependent on external resources and communications. During a volcanic crisis, the Volcanological and Seismological Observatory of Martinique (OVSM), located 8.5 km away from Mount Pelée, is crucial as its main goals include real-time monitoring of Mount Pelée and regional seismicity. The OVSM also plays a key part in crisis management by advising the competent authorities. Two additional major infrastructures include the Bellefontaine powerplant (located at 16.5 km from the volcano summit) that produces ≈ 60% of the total electric power in Martinique, and the Aimé Césaire international airport. The latter is located close to Fort-de-France (30 km from Mount Pelée) and would be essential both for the emergency response (human reinforcements and supplies sent from Guadeloupe and/or from the mainland) and for evacuation procedures. A detailed risk assessment is beyond the scope of this work, but we propose here to use our results to discuss the probabilities of critical ash thresholds being reached at these key locations.

Figure 10 shows the predicted area covered by ash as a function of threshold deposit thickness, determined from the iso-probabilities of 5, 10 and 50% shown in Figs. 8 and 9. The purple arrows and associated dotted lines show the circular area whose radius corresponds to the distance between each key location and the Mount Pelée summit. These areas and the associated probabilities of reaching a given deposit thickness are summarized in Table 3. Our calculations show that the OVSM could be impacted by ash thicknesses ranging from 2.5 cm (50% probability) to 84 cm (5% probability), with important consequences ranging from possible water contamination and extensive repair needed on supply networks to possible roof collapse and burial of some instruments. The Bellefontaine powerplant, further south, could be covered by 4 mm (50% probability) to 13 cm (5% probability) of ash, with associated damages ranging from ash infiltration and possible abrasion of some moving parts of gas turbines, to destruction of exposed equipments. In any case, such thicknesses could lead to temporary service disruption (Wilson et al. 2017). Finally, the international airport of Fort-de-France could be covered by an ash thickness of less than 1 mm (50% probability) to 2.5 cm (5% probability). There is thus a 50% probability of no disruption at all at the airport, but a non-negligible 5% probability of an airport closure and of possible abrasion of runway after cleaning operations (Wilson et al. 2017). It is also important to note that an ash thickness greater than 5 cm leads to extensive damage to most components of heating, air conditioning and ventilation systems (Wilson et al. 2017), which are often of paramount importance for these three key infrastructures.

Fig. 10

Predicted area covered by ash as a function of threshold deposit thickness, determined from Figures 8 and 9. The 5% (diamonds), 10% (squares) and 50% (triangles) probabilities of reaching the 1 mm (1.07 kgm⁻²), 1 cm (10.7 kgm⁻²), 10 cm (107 kgm⁻²), 20 cm (214 kgm⁻²) and 1 m (1,070 kgm⁻²) static thresholds are shown, as well as their power-law best-fit curve (black solid lines). The purple arrows and dotted lines show the circular area whose radius corresponds to the distance between Mt. Pelée and key locations in Martinique. The corresponding probabilities of reaching a given threshold thickness are given in Table 3 and discussed in the main text

Table 3 Probabilities of reaching critical thresholds at key locations in Martinique, according to Fig. 10

In order to investigate more distal consequences of Plinian eruptions, we consider two locations further south: the towns of Sainte-Anne in Martinique and Castries in Saint Lucia. The south of Martinique, represented by Sainte-Anne, could thus expect less than 1 mm (50% probability) to 4 mm (5% probability) of ash in case of a future Plinian eruption of the Mount Pelée. These thicknesses correspond to a range of damage from minor discomfort to required maintenance on supply networks. Note that during the phreatic eruption of May 2, 1902, 1 mm of ash was measured in Le Marin (close to Sainte-Anne) (Lacroix 1904), an observation that is in good agreement with our predictions. As for Castries in Saint Lucia, less than 1 mm to 1 mm of ash could be deposited in case of an eruption of Mount Pelée. Even if such thicknesses correspond to minor damage, it is a reminder that the consequences of explosive eruptions are of regional scale, hence that inter-islands cooperation and communication are necessary to ensure the safety of populations. We finally note that the area corresponding to the Fort-de-France international airport in Table 3 equates to the area encompassing the south of Dominica. The possible damages described in the previous section thus apply to this island as well. These results must however be taken with caution since the probability of northerly winds over Martinique can be significantly increased when a hurricane is passing by Martinique (Michaud-Dubuy et al. 2019). Such winds could spread ash over long distances, far beyond the south of Martinique. This result illustrates that multi-hazards need to be incorporated in a future integrated hazard map, along with other volcanic hazards described in the next section.

Other hazards

Although this work is dedicated to tephra fallout hazard assessment, six other volcanic hazards should be considered in Martinique, as suggested by Stieltjes and Mirgon (1998): pyroclastic flows, lava intrusions/flows, gas emissions, lahars, landslides and tsunamis (we could also add volcano-tectonic earthquakes induced by magma injection or withdrawal). In order to produce a new integrated volcanic hazard map for Mount Pelée volcano, a complete re-assessment of each hazard should be performed in the future. We argue that the Plinian eruptive history of Mount Pelée is much richer than previously thought, but many phreatic and Pelean eruptions certainly remain unknown as well. Since Mount Pelée exhibits some signals of increased activity and as hazard assessment strongly relies on eruption frequency and intensity, a careful revisit of this eruptive history is paramount.

An interesting example is the Pelean event named “Nuée de Balisier-Calave” (or NBC) by Traineau (1982). Our field investigations suggest that the main pyroclastic flow produced a secondary plume that impacted an area considered as safe in the current hazard map. This finding shows that even a Pelean eruption can have a strong impact far beyond the source if such a secondary plume forms. Such a phenomenon was often observed in similar volcanic context (e.g., the 1991 Unzen eruption, Watanabe et al. (1999); or the 2006 Tungurahua event, Engwell and Eychenne (2016)).

Conclusion

In this work, we proposed an integrated approach to refine tephra fallout hazard assessment in Martinique, based on eruptive scenarios determined from our revisited Plinian eruptive history of Mount Pelée volcano, and considering daily wind variability for the first time.

We first reproduced the VEI5 79 cal CE P3 eruption, the maximum expected eruption from our stratigraphic record, with the HAZMAP model. We then built a first hazard map for tephra fallout based on this scenario, which shows that even the south of Martinique would be impacted by a future powerful eruption of Mount Pelée volcano. We then produced a matrix of 16 scenarios based on the geological records and performed 239,600 simulations to elaborate the final hazard maps. We showed that the hazard for tephra fallout strongly depends on the eruptive scenario chosen and that the northern part of Martinique is strongly threatened regardless of the eruptive scenario. This conclusion is consistent with the current hazard map. However, the southern part of the island has a non-negligible probability of being impacted by a future eruption, both by tephra fallout and by its high dependance to the water and electrical network carried from the northern part of Martinique. Finally, we determined the range of tephra thicknesses that could be expected at key infrastructures in Martinique: the Volcanological and Seismological Observatory of Martinique (OVSM), the Bellefontaine powerplant, and the international airport, and showed that they have a non-negligible probability to be (sometimes strongly) impacted by a future Plinian eruption of the Mount Pelée. These results are consistent with our field measurements. Yet, both the southern part of Martinique and the airport are in areas considered as safe in the current hazard map.

These new hazard maps focus only on tephra fallout. To go even further, it is necessary to revisit the phreatic and Pelean eruptive history of Mount Pelée to re-assess the corresponding hazards and move towards a new integrated volcanic hazard map in Martinique. It is important to remember that these hazard maps are not risk maps, which would require an additional vulnerability assessment of the elements that may be affected during an eruption (e.g. population, buildings, networks), but can already be considered as a useful tool as Mount Pelée volcano has now entered into a new phase of activity. The new hazard maps presented in Figs. 5 and 9 will be included in their current form in the revised version of the French emergency plan (ORSEC) in response to a Plinian eruption in Martinique. This ongoing work performed in collaboration with the Prefecture of Martinique and DEAL Martinique will serve as a basis for the risk management and evacuation plan update in case of a future volcanic disaster.