Long-term hazard assessment of explosive eruptions at Jan Mayen (Norway) and implications for air traffic in the North Atlantic

Volcanic eruptions are among the most jeopardizing natural events due to their potential impacts on life, assets, and the environment. In particular, atmospheric dispersal of volcanic tephra and aerosols during explosive eruptions poses a serious threat to life and has significant consequences for infrastructures and global aviation safety. The volcanic island of Jan Mayen, located in the North Atlantic under trans-continental air traffic routes, is considered the northernmost active volcanic area in the world with at least five eruptive periods recorded during the last 200 years. However, quantitative hazard assessments on the possible consequences for the air traffic of a future ash-forming eruption at Jan Mayen are nonexistent. This study presents the first comprehensive long-term volcanic hazard assessment for the volcanic island of Jan Mayen in terms of ash dispersal and concentration at different flight levels. In order to delve into the characterization and modeling of that potential impact, a probabilistic approach based on merging a large number of numerical simulations is adopted, varying the volcano's eruption source parameters (ESPs) and meteorological scenario. Each ESP value is randomly sampled following a continuous probability density function (PDF) based on the Jan Mayen geological record. Over 20 years of meteorological data is considered in order to explore the natural variability associated with weather conditions and is used to run thousands of simulations of the ash dispersal model FALL3D on a 2 km resolution grid. The simulated scenarios are combined to produce probability maps of airborne ash concentration, arrival time, and persistence of unfavorable conditions at flight levels 50 and 250 (FL050 and FL250). The resulting maps can serve as an aid during the development of civil protection strategies, to decision-makers and aviation stakeholders, in assessing and preventing the potential impact of a future ash-rich eruption at Jan Mayen.

1 Introduction

Along with earthquakes, tsunamis, and weather extremes, explosive volcanic activity is among the most threatening natural hazards, with the potential to contribute to global warming and environmental changes (Ward, 2015). The impacts of volcanic emissions can extend over large distances from the source, posing a threat to human health and jeopardizing air navigation. Some recent examples of events leading to losses worth millions of US dollars due to air traffic disruption include the eruptions in Eyjafjallajökull (Iceland, 2010), Grímsvötn (Iceland, 2011), and Puyehue-Cordón Caulle (Chile, 2011) (Mazzocchi et al., 2010; Oxford-economics, 2010; Tesche et al., 2012; Karlsdóttir et al., 2012; Budd et al., 2011; Elissondo et al., 2016). These events were a stark reminder of the importance of volcanic hazard assessment and the related quantification of impacts of future eruptions, both essential tools to inform governments, aviation stakeholders, and society in general, contributing, in this way, to their preparedness. However, inferring what will be the impact of the next eruption is still a challenge for volcanologists, especially in the case of poorly characterized volcanoes. Classically, this can be done in two different ways: (1) by using the available information on the eruptive history of the volcano or (2) by using monitoring data and direct observations of the ongoing phenomena during a volcanic crises. The first case is what we call long-term hazard assessment, usually intended for cost–benefit analysis, long-term planning, and mitigation action design (Marzocchi et al., 2006), whereas the second one is a short-term hazard assessment that eventually would become a more deterministic forecasting tool for supporting the definition of emergency procedures.

In 2019, before the COVID-19 pandemic aviation break, Icelandic airports received around 8 million passengers (7 million international and 0.7 million domestic) on a total of 181 000 flights (Isavia, 2019). In turn, polar air traffic routes had shown a marked increase over the last few years, with a 15-fold increase between 2003 and 2015, reaching more than 14 000 flights a year from 2016 (NavCanada, 2017; Stewart-Green, 2016).

Although Jan Mayen (JM) volcano tephrochronology reveals at least eight eruptive periods over the last 600 years, five of them concentrated in the last 200 years (Gjerløw et al., 2016), the potential impact on air traffic following a future ash-forming eruption has never been assessed. According to Gjerløw et al. (2016), the most likely volcanism at JM is characterized by effusive Hawaiian to violent Strombolian eruptions and, to a lesser extent, by lava domes and Surtseyan eruptions. However, due to the possibility of magma interacting with seawater, snow, or ice, the likelihood of moderately to highly explosive eruptions is considerable. Historical distal records of trachytic tephra found in Ireland (Hunt, 2004) and basaltic tephra found in older sedimentary records in the North Atlantic (Lacasse and Garbe-Schönberg, 2001; Brendryen et al., 2010; Voelker and Haflidason, 2015) and in Greenland ice cores (Abbott and Davies, 2012) show the potential for producing Plinian explosive eruptions, whose size and frequency are, however, highly uncertain.

This paper presents the first comprehensive long-term probabilistic volcanic hazard assessment (PVHA) for the volcanic island JM focused on the potential impact of the airborne tephra concentration on Arctic and North Atlantic air routes. This is done by using the FALL3D model (Folch et al., 2009, 2020) to simulate the transport of ash clouds and their concentration at two relevant flight levels over a geographical area of approx 2000 km × 2000 km covering Iceland and the UK.

To account for the natural variability in volcanic eruption intensity, vent position, and wind field, two main steps have been followed as suggested in Sandri et al. (2016).

Firstly, based on field data, the possible eruptive classes for the JM volcano are identified and the probability density function (PDF) describing the relative probability of the different classes occurring is defined. For each class, a PDF for each eruptive parameter (such as eruption duration or total erupted mass) is defined in order to account for the natural variability in the eruption conditions. Then, by randomly sampling these PDFs, a large dataset of eruption source parameters (ESPs) to be used as model input is generated. A novel strategy is proposed to treat and describe the styles of pulsating eruptions, characterized by a series of discrete short-lived events followed by occasional interruption of tephra emission. Secondly, to fully explore the natural variability in the meteorological conditions, the numerical simulations have been randomly initialized within the period 1999–2020 (21 years). The meteorological data have been obtained from the ERA5 reanalysis dataset and FALL3D used to generate thousands of simulations per representative eruptive scenario. As a result, the following questions are answered:

• What is the probability that, in the case of an eruption at JM, the ash cloud concentration will exceed the critical conditions for safe flights within a domain extending down to the UK airspace at 3, 6, 12, 18, and 24 h after the beginning of the eruption?

• In the case of an eruption at JM, what is the probability that airports in Iceland and the UK will be affected by the presence of ash?

• What is the probability of exceeding a predefined hazardous temporal persistence of unsafe flight conditions?

• Which flight level (FL) is likely to be predominantly affected by critical concentrations of volcanic ash?

The rest of the paper is organized as follows: Sect. 2 provides a historical overview of the Holocene volcanic activity of the JM volcano. Section 3 describes the most likely eruptive classes based on the five known historical eruptions of JM, fits them into a PDF for the total erupted volume, addresses a novel strategy to treat them, and describes the styles of pulsating eruptions. Sections 4 and 5 present results and discussions. Finally, Sect. 6 concludes the study.

2 Jan Mayen volcanism

JM is a Norwegian volcanic island located in the North Atlantic Ocean at 71^∘ N, 8^∘ W around 600 km north of Iceland, in the Norwegian Greenland Sea (Fig. 1). According to Kandilarov et al. (2012) and Larsen et al. (2021), the Jan Mayen microcontinent (JMMC, Fig. 2a) is a structural entity enclosing the Jan Mayen Ridge (JMR) and the surrounding area, including the Jan Mayen Basin (JMB), the Jan Mayen Basin South (JMBS), the Jan Mayen Trough (JMT), and the Southern Ridge Complex (SRC) (see Fig. 2b). To the north, the JMMC is bordered by the Jan Mayen Fracture Zone (JMFZ) and the volcanic complex of Jan Mayen, while to the south, east, and west it is bordered by the NE coastal shelf of Iceland (NIS), the Norway Basin, and the Kolbeinsey Ridge (KBR) respectively (Fig. 2a). Although the historical activity reports at least five eruptive periods over the last 200 years (since the discovery of the island at the beginning of the 17th century) (Gjerløw et al., 2016), its Holocene eruptive history is basically unknown. In this sense, the eruptive history of JM comprises only a very few distal sediment cores as well as lava flows and tephra deposits from eruptions on the ice-free parts of the Beerenberg volcano. Distal records of trachytic tephra found in Ireland (Hunt, 2004) and basaltic tephra found in older sediment records in the North Atlantic (Lacasse and Garbe-Schönberg, 2001; Brendryen et al., 2010; Voelker and Haflidason, 2015) and in Greenland ice cores (Abbott and Davies, 2012) have shown the potential for explosive ash-forming eruptions, whose size, frequency, and potential impact are, however, uncertain. According to Imsland (1978), explosive hydromagmatic eruptions were common earlier in the history of JM. Nevertheless, as the island grew above sea level, such eruptions became less frequent and the volcanism essentially localized in two different regions: (1) the Beerenberg central volcano and its flank eruptions in the northeastern part (called Nord-Jan) and (2) the volcanic ridge extending from the middle to the southwestern part (called Midt-Jan and Sør-Jan respectively). On one hand, considering that the higher altitudes of the volcano are ice-covered and glacier tongues extend down to sea level at several locations, the Holocene eruptions from the summit crater are difficult to map, and no tephra layers have been positively linked to eruptions from the summit. Only a few land-based tephra records on the ice-free areas of Beerenberg have been mapped with some detail. Based on several sediment cores, Gjerløw et al. (2016) conclude that the Holocene volcanism on Beerenberg has been effusive or mildly explosive. As a result, the most common forms of recognized volcanic activity at Beerenberg are flank eruptions in the form of basaltic fissure and Strombolian to violent Strombolian eruptions. Eruption frequency is difficult to assess due to scarce reconstruction data. However, during historical times, the Beerenberg's eruption rate was around one eruption every 60–70 years, with eruptive phases lasting in the range of days to months. During the most recent effusive eruption in 1970, the largest known one during the Holocene, the volume of lava flows was at least 0.5 km³ dense rock equivalent (DRE) (Siggerud, 1972). On the other hand, volcanism on Midt-Jan and Sør-Jan represents mostly effusive eruptions characterized by scoria cones, shallow marine to coastal phreatomagmatic eruptions, coulees, and domes (Larsen and Guðmundsson, 2016; Gjerløw et al., 2015). The eruption frequency on this part of JM is also difficult to assess due to erosion and the superimposition of newer vents (possibly covering and removing older ones). However, considering visible evidence, the (under) estimated number of eruptions over the last 10 000 years is around 45, resulting in an eruption frequency of 1 eruption every 220 years. The duration of the eruptions from Midt-Jan and Sør-Jan is still unknown. The unrest episode recorded in 1732 (Eggøya, Midt-Jan), which led to the largest known explosive eruption, was a Surtseyan eruption that dispersed tephra over large parts of JM and the surrounding seas. The volume of tephra ranges between 0.3–0.4 km³ (Gjerløw et al., 2015).

Figure 1 JM location and computational domain for the JM PVHA including Iceland, Ireland, and the United Kingdom (blue box). The red contour shows the FIR (flight information region) for which the Icelandic Meteorological Office is responsible (for visualization purposes only). The blue star and triangle on the zoomed-in map indicate the location of the Beerenberg volcano and Eggøya crater (1732 Surtseyan eruption) respectively. The two blue circles show the two hypothetical vent locations in the wind profile analysis. Black circles correspond to Keflavik and Akureyri (Iceland), Vágar (Faroe Islands), Edinburgh (Scotland), and Heathrow (London, UK) airports. Nord-Jan, Midt-Jan, and Sør-Jan correspond to the three areas into which JM is classically divided.

Figure 2 Overview map (a) of the study area with the location of structural elements identified in potential field data. Structural elements map (b) for the Jan Mayen microcontinent (JMMC): mapped faults, fractures zones, and lineaments based on Peron-Pinvidic et al. (2012) and Gernigon et al. (2015). The background image is shaded bathymetry (IBCAO 3.0; Jakobsson et al., 2012; Amante and Eakins, 2009). Image retrieved from Blischke et al. (2017).

3 Methodology

3.1 Eruption scenarios

Despite the limitations of a complete geological record composed of both chronological and statistical data, the possible relative eruptive scenarios at JM are based on five known historical and prehistorical eruptions. According to the categorization proposed by Larsen and Guðmundsson (2016) and Gjerløw et al. (2016), eruption scenarios can be characterized by small (<0.1 km³), moderate (0.1–0.5 km³), and large (>0.5 km³) DRE volumes or magnitudes (see Table 1) (Pyle, 2015) and sub-Plinian eruptions.

• Small eruptions are mostly effusive events characterized by small lava flows or small scoria cones, with erupted volumes ranging 10⁷–10⁸ m³ (total less than 0.1 km³ DRE), corresponding to eruption magnitudes of 1 to 2; hence the volcanic explosivity index (VEI) is 2 (Newhall and Self, 1982). Based on historical occurrences, this scenario can last for about 35–40 h.

• Medium eruptions include subaerial, subglacial, and even Surtseyan eruptions depending on within which environment they occur. Subaerial eruptions would be mainly located on Beerenberg volcano, and they are expected to be effusive and/or Vulcanian to violent Strombolian. When effusive, medium eruptions are characterized by not only aa lavas but also pahoehoe flows. Surtseyan eruptions are expected to be located on JM and in the surrounding shallow part of the ocean. These eruptions consist of phreatomagmatic pulses, each of which, according to observations, can last for approximately 0.5–8 d, generate volcanic plumes between 3 and 11 km a.s.l. (above sea level), and have a range of total erupted volume of 10⁸–10^8.7 m³ (0.1–0.5 km³ DRE), corresponding to eruption magnitudes 3 to 4, and VEI = 3. The total duration of the eruption is not well constrained, as it can last between approximately 4 d and 1 month. As a result, tephra-forming phases are expected, producing deposits more than 1 m thick within 5 km of the vent. The reference eruption for the Surtseyan type is the 1732 CE Eggøya eruption that produced at least 0.3–0.4 km³ of tephra (0.16–0.21 km³ DRE) (Gjerløw et al., 2015, 2016).

• Large eruptions are expected to be initially subglacial and include moderate to sub-Plinian eruptions. During the opening phases, due to magma–ice interaction, the activity is explosive and characterized by plume heights reaching more than 10 km a.s.l. and a range of total erupted volume of 10^8.7–10⁹ m³ (total volume emitted > 0.5 km³ DRE), corresponding to eruption magnitudes of 4 to 5, and VEI = 4. In this initial short-lasting explosive phase, a very small amount of tephra is expected to be ejected (approximately 5 % of the total erupted mass). The reference eruption for this type is the 1970 event that produced at least 0.5 km³ DRE (Siggerud, 1972). As the eruption proceeds, it becomes more effusive, lasting for 1–4 d.

• Very large eruptions include sub-Plinian to Plinian eruptions characterized by column heights from 15 km to 25 km a.s.l. and a range of total erupted volume of 10⁹–10^9.7 m³, hence with eruption magnitudes of 5 to 6, corresponding to sub-Plinian type I or VEI ≥ 5. According to Gjerløw et al. (2016), there is evidence in geological records (extending beyond the Holocene) of 10 tephra layers from sub-Plinian and/or Plinian events in 119 000 years. Because of this, we assign a subjective probability of 1 % to this class in the case of an eruption.

Table 1 Possible relative eruption scenarios on JM. The categorization is based on the volume of tephra emitted in DRE. Data obtained from Larsen and Guðmundsson (2016) and Gjerløw et al. (2016). According to the geological record (extending beyond the Holocene), sub-Plinian/Plinian events are highly unlikely (1 %). Because of this, they are not included in this table.

3.2 Probabilistic hazard assessment approach

Until a few years ago, volcanic hazard assessment was largely based on the concept of the “eruptive scenario”, characterized by subjectively defined eruption conditions. Hazard was then quantified under the strong assumption that the next eruption from a given volcano would be similar to the selected “representative eruptive scenario” (Macedonio et al., 2008; Barsotti et al., 2018). However, when assuming a representative eruptive scenario, one is implicitly neglecting the large uncertainties (both aleatory and epistemic) in the parameters that define the scenario, also called “intra-size variability” (Woodhouse et al., 2015; Harvey et al., 2018). More recent approaches try to circumvent the effects of natural variability by averaging hundreds of simulations where eruption parameters are sampled within a broad set of eruptive conditions in the so-called “eruption range scenarios” (Bonadonna et al., 2005; Folch and Sulpizio, 2010; Prata et al., 2019). However, the use of a specific and limited range of eruption parameters continues to introduce a large bias and uncertainties into the description of potential eruptive processes. For this reason, more recent approaches are based on the concept of a continuum of possible combinations of eruptive parameters, which translates into exploring a large set (many thousands) of simulations as proposed by Sulpizio et al. (2012) and Sandri et al. (2016). Eruption parameters (e.g., total erupted mass, duration of the fallout phase, mass eruption rate, total grain size distribution) are defined and randomly sampled from specific probability distributions (Sandri et al., 2016). The processes for sampling and weighting possible statistical combinations of values for the volcanological parameters correspond to their probability of occurrence: this allows giving more/less weight to more/less likely combinations. In order to explore the intra-size variability, we proceed as in Sandri et al. (2016):

1. A very broad range of possible eruptive scenarios, characterized by the total erupted volume, are selected as explained in Sect. 3.1.

2. The total erupted volume is used to define the total erupted mass, the eruption magnitude, and the VEI.

3. The eruptive range is split into eruption classes linked to representative members (see Sect. 3.1), each characterized by an approximate conditional probability in the geological and historical record (see Table 1).

4. Over the total range of possible erupted volumes (approximately 10⁷–10¹⁰ m³), up to six different truncated probability density functions (PDFs) are tested to describe the conditional probability of these four mutually exclusive classes: normal, exponential, log-logistic, lognormal, gamma, and Weibull. The best model is selected according to Akaike's information criterion (AIC) (Bozdogan, 1987; Akaike, 1998), where the relative goodness of fit of such PDFs (i.e., the likelihood of the catalog's frequencies of the different eruptive classes, under different PDFs) is compared, penalized by the number of parameters. The PDF with the lowest AIC is considered the best among all the specified ones. Indeed, the assumption of a common PDF for the total erupted volume across the different eruption classes allows smooth and coherent linking among them (Sandri et al., 2016). For JM, the Weibull PDF better fitted the expected frequencies on the sub-ranges for the four different eruption classes. This PDF has been used to assign a conditional probability of occurrence to each simulation as a function of the associated total erupted volume (see Fig. 3).

5. Considering the behavior of similar scenarios including wet plumes, for medium and large classes we account for ash aggregation assuming two different aggregate types characterized by densities in the range of 250 and 350 kg m⁻³ and diameters between 100 and 250 µm.

Figure 3 Weibull PDF describing the conditional probability of different eruptive magnitudes in the case of an eruption at JM. The four colors cover the erupted volume ranges in the four “classical” eruption classes for JM, classically synthesized into four representative scenarios with a fixed mass, neglecting the variability in volume around these scenarios. The area under the different parts of the plot corresponds to the probability of an effusive, medium, large, and sub-Plinian class range eruption respectively, conditional to eruption occurrence. These values are in agreement with previous studies for JM (Larsen et al., 2017; Gjerløw et al., 2016).

Figure 4 Proposed strategy to treat and describe the styles of pulsating eruptions, characterized by a series of discrete short-lived events followed by occasional interruption of the tephra emission.

3.3 Pulsating eruptions – modeling and strategy

A novel strategy is proposed to describe the styles and model dispersal from pulsating eruptions (Surtseyan eruptions, belonging to the medium class in Table 1), characterized by a series of discrete short-lived events followed by occasional interruption of the emission of tephra. The strategy has been developed considering the ranges of all the ESPs described in Sect. 3.1. For each pulsating scenario, the ESPs associated with the column shape, total grain size distribution, and sphericity of tephra particles are also sampled from given PDFs. However, the difference is that column heights are not derived from the mass eruption rate but using the following approach (see Fig. 4).

1. Random sampling is undertaken of both the total erupted mass (TEM_T) and the total duration of the eruption (Dur_T) considering values reflected in Sect. 3.1.

2. If the sum of masses erupted by all pulses does not equal or exceed the total erupted mass previously sampled, loop to the following:

   • Create the ith pulse, sampling randomly column height (H_i) and duration (Dur_i). The duration is sampled from a normal distribution consistent with the data reflected in Table 1. The column height is sampled from a triangular distribution with the lower limit at 3 km, peak at 6 km, and upper limit at 11 km a.s.l.

   • Compute the total erupted mass for such a pulse (TEM_i) using the Mastin et al. (2009) relationship.

   • Compute ($\sum _{i=\mathrm{1}}^n$ (TEM_i)), with n being the number of pulses generated so far.

     • If $\left(\sum _{i=\mathrm{1}}^n\left({\mathrm{TEM}}_i\right)\right)$ > 0.97 ⋅ TEM_T and $\left(\sum _{i=\mathrm{1}}^n\left({\mathrm{TEM}}_i\right)\right)$ < TEM_T, modify TEM_i to obtain $\sum _{i=\mathrm{1}}^n$ (TEM_i) = TEM_T, thereby avoiding small pulses. Compute the new column height (H_i) using Mastin et al. (2009).

     • Else, if total erupted mass obtained $\left(\sum _{i=\mathrm{1}}^n\left({\mathrm{TEM}}_i\right)\right)$ < TEM_T, save the pulse. Otherwise, discard the pulse.

3. Compute the duration of the eruption as the sum of the duration of all the pulses $\left(\sum _{i=\mathrm{1}}^n\left({\mathrm{Dur}}_i\right)\right)$. If $\sum _{i=\mathrm{1}}^n\left({\mathrm{Dur}}_i\right)<{\mathrm{Dur}}_{\mathrm{T}}$, generate n−1 inter-pulses at a random time (Res_i) so that their sum equals δ ($\mathit{\delta }={\mathrm{Dur}}_{\mathrm{T}}-\sum _{i=\mathrm{1}}^n\left({\mathrm{Dur}}_i\right)$) and insert them between pulses. Otherwise, if $\sum _{i=\mathrm{1}}^n\left({\mathrm{Dur}}_i\right)>{\mathrm{Dur}}_{\mathrm{T}}$, update Dur_T to equal $\sum _{i=\mathrm{1}}^n\left({\mathrm{Dur}}_i\right)$. This case actually supposes a continuous eruption where each pulse occurs without a rest period.

3.4 Vent location sensitivity

Given the scales of JM and the considered domain, the effects of the uncertain vent location on the resulting long-range hazard assessment can be expected to be negligible. In order to check this assumption, we inspected how ERA5 wind profiles vary along the island by focusing on two vent locations at the NE (71.15^∘ N, 7.95^∘ W) and SW (70.82^∘ N, 9.02^∘ W) edges of the island, approximately 55 km apart (blue circles in Fig. 1 inset). At these locations we inspected the following:

• Local wind profiles. Figure 5 shows vertical profiles of the wind speed and direction averaged for the whole month of December 2019. As observed, there are few differences in patterns between the two locations.

• Annual wind profiles. Figure 6 shows the wind profiles averaged monthly for the year 2018. Once again, there are no differences between the two locations.

Figure 5 Monthly averaged ERA5 wind profiles (speed and direction) at two different locations (NE and SW) of JM.

Figure 6 Monthly comparison of the wind pattern computed in two different locations (NE and SW) for a whole year, 2018. (a, b) Wind patterns corresponding to the NE vent. (c, d) Wind patterns corresponding to the SW vent. Results were obtained by averaging 1 year of ERA5 data.

Considering the current limitation of both the grid resolution and the meteorological data resolution, the location of potential JM vents does not influence the ash dispersal pattern. As a result, we will not consider the uncertainty of the vent location and we assume a fixed vent in the middle of the island.

4 Results – hazard maps and uncertainty quantification

Hazard and probability maps (Elefante et al., 2010) are powerful tools to inform users about the spatial and temporal potential impact of specific volcanic phenomena. Commonly, they consist of exceedance probability curves, referred to as hazard curves (Hill et al., 2013). These hazard curves quantify, in a grid point of the target domain and within a specific time window (exposure time) (Budnitz et al., 1997), the exceedance probability of an intensity measure threshold for a specific phenomenon (e.g., tephra load at ground level or airborne tephra concentration).

Our objective is to show the usefulness of probabilistic volcanic hazard assessment in the framework of high-performance computing, evaluating the impact of low-probability but high-consequence events on air traffic (between Iceland and the UK; see Fig. 1) from a potential eruption at JM while quantifying how the ESPs and wind patterns (velocity and direction) influence hazard and probability maps of ash dispersal and airborne tephra concentration.

Figure 7 Overview of volcanic ash or sand/dust impacts on jet engines. Modified from a Rolls-Royce review of aircraft encountering airborne particle clouds (Ellis et al., 2021; Clarkson et al., 2016).

Although our method allows analyzing any desired FL, in this work, only FL050 and FL250 (approx. 1.5 and 7.5 km a.s.l.) will be analyzed. Such FLs were selected considering standard cruise and maximum risk action altitudes like takeoff or landing. Finally, three selected ash concentration thresholds (0.2, 2, and 4 mg m⁻³) were selected based on the impacts of volcanic ash on jet engines summarized in Fig. 7 and the considerations included in the Volcanic Ash Contingency Plan published by the International Civil Aviation Organization (ICAO, 2021). As a result, we analyze the results using isolines at FL050 and FL250 for these three selected ash concentration thresholds through three different types of probabilistic map:

• arrival time maps – the expected time required for the ash concentration to exceed 2 mg m⁻³ at FL050 with an exceedance probability of 5 %, between 0 and 48 h since the beginning of the eruption;

• exceedance probability maps – reporting the probability of reaching ash concentration above a given threshold (0.2, 2, and 4 mg m⁻³) at FL050 and FL250 at some time from the onset of the eruption to up to 48 h after its end;

• persistence maps – showing the fraction of hours (since the beginning of the eruption) during which the ash concentration exceeds a given threshold (0.2, 2, and 4 mg m⁻³) with a probability larger than 5 %.

Figure 8 The isolines show the arrival time in hours required for the ash concentration (at FL050) to exceed a threshold of 2 mg m⁻³ with an exceedance probability of 5 % between 0 and 48 h after the eruption. Black circles correspond to Keflavik and Akureyri (Iceland), Vágar (Faroe Islands), Edinburgh (Scotland), and Heathrow (London, UK) airports.

Figure 8 depicts the arrival time maps for large and medium eruptions respectively. The percent value in exceedance probability has been subjectively selected. However we highlight that our method allows a potential end-user to explore any value of exceedance probability: here, we only show the 5 % maps as an example.

Figures 9 and D1 in Appendix D show the probability of reaching or exceeding ash concentration above 0.2, 2, and 4 mg m⁻³ at FL050 and FL250 at some time from the onset of the eruption to up to 48 h after its end for the large and medium eruptive classes respectively.

Figure 9 Exceedance probability maps at FL050: the isolines show the probability of reaching an ash concentration above 0.2 mg m⁻³ (a, d), 2 mg m⁻³ (b, e), and 4 mg m⁻³ (c, f) at FL050 at some time from the onset of the eruption to up to 48 h after its end. The three upper plots apply to the large class and the three at the bottom to the medium one.

Thousands of eruptive scenarios have been simulated to reproduce and capture, in a probabilistic way, the variability in phenomena which can vary strongly in space and time. However, proper uncertainty quantification (UQ) is needed to quantify how reliable the prediction is. Such UQ can provide useful knowledge about the diversity of the dominant winds, the range in the airborne tephra concentration and its extent depending on the type of eruption, the ESPs related to the eruption size, and the feature of pulsating events for medium-size eruptions. As a consequence, the threat evaluation and the spatio-temporal analysis presented here could bring forth a more robust comprehensive hazard assessment.

Kristiansen et al. (2012) concluded that the main source of epistemic/aleatory uncertainty in ash dispersal forecasts comes from the quantification of the eruption source term (eruption column height and emission rate). Here, we address the quantification of uncertainty over the airborne tephra concentration and its extent. To do that, we assess the 95 % confidence interval (i.e., range between the percentiles of 97.5 and 2.5) in the probability distributions describing the hazard curves for the concentration of tephra for each point in the domain. These probability distributions are deeply related to the number of simulations or scenarios used that model such concentrations; a detailed analysis of how the number of simulations affects the sensitivity of this uncertainty can be found in Appendix A.

Figure 10 shows different maps, at different levels of confidence, produced by cutting the point-wise hazard curves at different percentiles. They correspond to a 4D analysis where concentrations are the highest at any time from the onset of the eruption to up to 48 h after its end.

Figure 10 Concentration extent hazard map at FL050: relative uncertainties related to airborne ash cloud concentrations above 0.2, 2, and 4 mg m⁻³ and extent. Each map corresponds to a different level of confidence, produced by cutting the point-wise hazard curves at different percentiles.

Figure 11 shows, from top left to bottom right, the probability of reaching or exceeding ash concentration above 2 mg m⁻³ at FL050 for more than 1, 3, 6, 12, 18, and 24 h respectively from the onset of a large eruption to up to 48 h after its end. Figure 12 shows the same but for the medium-size eruptions, and, in Appendix D, Figs. D2 and D3 display the same information as Figs. 11 and 12 but for FL250.

Figure 11 Persistence maps at FL050 (large class): the isolines show the probability of reaching or exceeding an ash concentration above 2 mg m⁻³ at FL050 for 1, 3, 6, 12, 18, and 24 h during the eruption to up to 48 h after its end.

Figure 12 Persistence maps at FL050 (medium class): the isolines show the probability of reaching or exceeding an ash concentration above 2 mg m⁻³ at FL050 for 1, 3, 6, 12, 18, and 24 h during the eruption and to up to 48 h after its end.

5 Discussion

Results have been carried out considering the intrinsic limitations of the methodology (partially due to the scarcity of data related to complete geological records composed of both chronological and statistical data). This fact is important, as future advances in the geological catalog could have implications for future work assessing volcanic hazards and mitigating measures on JM.

5.1 Arrival time maps

As shown in Fig. 8, an ash-rich eruption originating from the JM volcano has the potential to affect the air traffic over Iceland after 36 h and, to some extent, the Faroe Islands or archipelago 48 h after the beginning of the eruption, with an exceedance probability of 5 %.

Table 2 Airport locations, azimuth, and distance to JM. Locations are expressed in lat–long coordinates. Azimuth and distance correspond to the azimuth in degrees and the distance between JM and the different airports.

Figure 13 shows the evolution of the probability over time of ash concentration exceeding 2 mg m⁻³ at FL050 and the international airports of Keflavik and Akureyri (Iceland), Vágar (Faroe Islands), and Heathrow (London, UK) (see Table 2 for distance references). The probability at any airport is neglectable during the first hours (approx. 10 and 15 h for the medium and large classes respectively) and then increases until stabilizing after several days (approx. 7 and 5 d for the medium and large classes respectively). For both eruption classes, Vágar Airport has a higher probability of exceeding such a threshold than the other nearest airports such as Keflavik. This is due to a very marked difference in the wind patterns between the north-northeast and the west. At 48 h after the beginning of the eruption, only Akureyri Airport should exceed probabilities above 5 % to reach the concentration threshold of 2 mg m⁻³ for a medium-class eruption. No airport shows exceedance probabilities for such a critical threshold above 25 %.

Figure 13 Exceedance probabilities vs. arrival times required for the ash concentration (at FL050) to exceed a threshold of 2 mg m⁻³ at different airports for medium and large eruptions from the beginning of the eruption.

5.2 Exceedance probability maps

Figures 9 and D1 in the Appendix D show substantially different results regarding ash concentration and extent for both eruption classes. Regarding concentrations, for the large class, values above 2 mg m⁻³ (reaching even 4 mg m⁻³, originally considered a no-fly zone) would affect part of the Icelandic FIR with exceedance probabilities between 10 % and 50 %. Instead, for the medium class, such concentrations would affect only low flight levels. Above FL250, moderate–higher probabilities are only found for polar routes. This is because the height of the eruptive column for medium-magnitude-class eruptions does not exceed 11 km (see Sect. 3.1). These results are in agreement with those shown in Fig. 10 (where maps, at different levels of confidence, were produced by cutting the hazard curves at distinct percentiles to depict the relative uncertainties related to airborne ash cloud concentrations and extent for both eruption classes at FL050).

Concerning the extent, at FL050, ash concentrations of up to 2 mg m⁻³ could reach almost the entire Icelandic FIR with probabilities between 10 % and 50 % for both eruption classes. This could threaten the vast majority of flights to and from northern routes. In contrast, at FL250, ash clouds would affect the northeast of Iceland for concentrations of up to 0.2 mg m⁻³ only when a large eruption occurs. Thus, we can conclude that for the medium class, only polar routes above FL250 would be threatened.

Finally, in a similar way to other types of analysis such as tephra ground load and probabilistic seismic hazard assessment, Fig. 14 provides a graphical representation of relative uncertainties related to airborne ash cloud concentrations above 0.2, 2, and 4 mg m⁻³ at FL050 at Keflavik Airport. This result, computed at each point of the target domain, could be eventually used as input for risk analysis such as for producing fragility curves and tolerance analysis and in general investigation of impact on infrastructure. In this view, it represents the most complete way to quantify hazard. Specifically, no dramatic differences are seen depending on the eruption size, and there is a non-negligible probability of overcoming the 2 mg m⁻³ threshold, even for low percentiles, given an eruption.

Figure 14 Concentration analysis: relative uncertainties related to airborne ash cloud concentrations above 0.2, 2, and 4 mg m⁻³ at FL050 at Keflavik Airport.

5.3 Persistence maps

According to Fig. 7, jet engines exposed to concentration conditions of up to 4 mg m⁻³ for more than 3 h would require inspection. This motivates the spatio-temporal analysis of persistent high-concentration scenarios. In terms of the ash cloud extent, results (Figs. 11 and 12) are slightly different: at FL050 a large part of the Icelandic FIR (reaching to some extent the Faroe Islands) with probabilities between 10 % and 50 % would be affected for both eruption classes for up to a total of 6 h after the beginning of the eruption. Instead, at FL250 (Figs. D2 and D3), such conditions would affect only high-latitude air routes (above 68^∘ N).

For high-concentration scenarios lasting longer than 12 h, some differences between the eruption classes can be observed. At FL050, an ash cloud has probabilities lower than 10 % of reaching latitudes as low as 68 and 66^∘ N for large and medium eruption respectively. These southernmost latitudes increase for higher persistence values, meaning that (obviously) only closer to the source may we get long-persisting clouds.

Finally, Fig. 15 presents a persistence analysis for the airports considered in this study, showing the exceedance probability of reaching ash concentration above the critical condition for maximum risk actions like takeoff or landing until 24 h after the beginning of the eruption. The most affected airports are Akureyri, Vágar, and Keflavik. London Heathrow (LHR) has very low probabilities (1.5 %–2 %) associated with 1 to 6 h of persistence scenarios until 48 h after the beginning of the eruption. Concerning the level of persistence, both eruption classes have similar behavior. Scenarios with persistence longer than 18 h are highly unlikely. However, when analyzing probabilities, the medium class reaches values twice as high as the large one. This observation can again be associated with their eruptive dynamic. The sustained injection of tephra into the atmosphere related to a series of discrete short-lived events increases the probability of prolonged high-concentration scenarios.

Figure 15 Persistence analysis for ash concentration of up to 2 mg m⁻³ at FL050. Exceedance probability of reaching or exceeding ash concentration above 2 mg m⁻³ at FL050 at 1, 3, 6, 12, 18, and 24 h at different airports during the eruption to up to 48 h after its end.

6 Conclusions

Despite the fact that the limitations of the geological records may constitute a bias for long-term volcanic hazard assessment, this work provides the first comprehensive analysis focused on the potential impact of airborne tephra concentration on Arctic and North Atlantic air routes due to an ash-forming eruption at Jan Mayen.

By looking at the field data, the possible eruptive classes for the JM volcano and their relative occurrence probability (expressed as probability density functions) have been identified. Here, medium to large eruptions (VEI of 3 to 4) are considered the most likely events, whereas Plinian eruptions are expected to occur with a probability of lower than 1 % (and for this reason they have not been considered in the analysis). A novel strategy to treat and describe the styles of pulsating eruptions (the case of the medium class), characterized by a series of discrete short-lived events followed by the occasional interruption of tephra emission, has been proposed. The natural variability in the eruption conditions for each possible eruptive class is described by different eruption source parameters. Such ESPs are randomly sampled from previously defined PDFs to generate a large dataset to be used as model input.

The results presented here are indicative of the potential impact of an eruption at JM in light of what is known about its previous eruptive activity. They could be used for land-use planning, defining mitigation actions, identification of vulnerable infrastructure, and cost–benefit analysis but not for emergency response.

Our results show that an ash-rich eruption originating from the JM volcano has the potential to affect air traffic over Iceland (after 36 h) and, to some extent, the Faroe Islands, after 48 h. Concerning airborne ash concentration and extent, for large eruptions, concentrations above 2 mg m⁻³ (even 4 mg m⁻³, originally considered a no-fly zone) would affect part of the Icelandic airspace (at different flight levels) with exceedance probabilities between 10 % and 50 % at some time from the onset of the eruption to up to 48 h after its end. For medium eruptions, these dangerous concentrations would affect only low flight levels (FL050). Above FL250 only polar routes would be affected.

When analyzing persisting concentration conditions where aircraft engines are exposed to a high concentration for more than 3 h, we conclude that at FL050 a large part of the Icelandic FIR (reaching to some extent the UK) would be affected for both medium and large eruption classes with probabilities between 5 % and 50 %. At FL250 such risky conditions would affect only high-latitude air routes (above 68^∘ N).

Finally, we want to highlight the robustness of our PVHA in terms of uncertainty quantification, which should be routinely considered in all these kinds of studies.

Appendix A: Simulation setup

This section provides an overview of the high-performance computing (HPC) environment used in this study and the setup process associated with the FALL3D model to simulate the eruptive scenarios. The most relevant settings to optimize the computational resources, as well as the simulation scheme followed, are described.

To account for the meteorological statistics in the simulation results of each eruptive class (large and medium), 1500 wind fields over the time period 1999–2020 were randomly sampled from the ECMWF ERA5 reanalysis dataset. Then, 1500 simulations (scenarios) combining meteorological conditions and ESPs for each class were run. Since the medium eruptive class is characterized by a series of discrete short-lived events, the total number of scenarios for this class was 3763. The total CPU–GPU hours used was 9.6 million, considering the architectures shown in Table C1 in Appendix C.

Figure A1 shows the variance of the tephra concentration at a given grid point with respect to the number of scenarios simulated. After 900 scenarios, the variance of the concentration begins to stabilize. This stabilization also suggests a reduction in uncertainty related to the intra-size variability in the eruptive scenarios themselves.

Figure A1 Value of the variance of tephra concentration for a given grid point with respect to the number of scenarios simulated.

All the scenarios were run using the FALL3D-8.1 model (Folch et al., 2020) over a 2000 km × 2000 km domain between 50 and 73^∘ N (latitude) and 2 and 24^∘ W (longitude) with a resolution of about 2 km. Eruptive vents were placed at 70.98^∘ N, 8.38^∘ W and 71.10^∘ N, 8.13^∘ W respectively for medium and large eruptive classes.

Appendix B: Sampling and processing workflow

Table B1 PDFs and values ranges of the main eruptive parameters for JM. Bounds on mass eruption rate values are a consequence of the sampling procedure for total erupted volume (Fig. 3) and duration of the fallout phase described in the text. For the total grain size distribution, references were chosen from the 1732 Eggøya and 2004 Grímsvötn eruptions for medium and large classes respectively. Erupted volumes are between 0.1–0.5 and >0.5 km³ for medium and large class ranges.

For each eruptive class, the PDFs describing each ESP were fixed following Sandri et al. (2016). It is important to note that this work only addressed the medium and large eruptive classes. Table B1 summarizes the PDFs and values ranges of the main ESP for JM. The sampling process can be described as follows:

1. Sample a value for the total erupted volume (or magnitude), duration of the fallout phase, column shape, and sphericity of tephra particles from their PDFs. The total erupted volume, expressed as DRE, is computed uniformly within a range of values (10⁸–10^8.9 m³). In a second step, we weight each total eruption volume based on the Weibull distribution function previously defined (Fig. 2). In this way, the unlikely events are properly represented.

2. Compute the mass fraction (%) associated with tephra fallout concerning the total erupted mass according to the available estimations from field data analysis. For the medium eruptive class, the value of the mass fraction is fixed to 0.8, whereas for the large eruptive class it is randomly sampled from [0.05, 0.10].

3. Compute the mass eruption rate and the column heights from the total erupted volume sampled. The mass eruption rates range between 3.009×10⁴–1.5×10⁶ and 6.94×10⁴–1.39×10⁶ kg s⁻¹ for the medium and large eruptive sizes respectively. The source term (vertical distribution of mass in the eruption column) is given by a Suzuki distribution (Suzuki, 1983) with the parameter A in the range of 3 to 4.5 and λ equal to 1.

4. Sample a time for the eruption start over a period of 21 years (1999–2020) considering the corresponding meteorological fields for the duration of the fallout phase and associate this randomly with a combination of the volcanological parameters. For this, ECMWF ERA5 reanalysis meteorological data have been used associated with the date and duration of the eruption. Modify FALL3D's input file with both meteorological and newly sampled data values.

5. Run FALL3D to obtain the tephra loading at different flight levels. The Ganser terminal velocity model (Ganser, 1993) and the CMAQ model parameterization for horizontal diffusion (Folch et al., 2009) were used as part of the model physics configuration.

6. Computet the outputs obtained from FALL3D. As a result, hazard and probabilistic maps describing the airborne ash concentration and time persistence at different flight levels on a large-scale and high-resolution domain are obtained.

Typical tephra particle densities and total grain size distributions were chosen to be consistent with previous values reported for the Eggøya (1790) and Grímsvötn (2004) eruptions for medium and large classes respectively. The type of aggregates was also chosen to be consistent with previous values reported for similar Surtseyan and phreatomagmatic eruptions.

Appendix C: Computational resource

Table C1 Joliot-Curie supercomputer. Characteristics corresponding with the two partitions available in this study.

Figure C1 Strong scalability analysis (time to solution). (a, b) Speedup and parallel efficiency analysis at Irene ROME (128 AMD processors per node). (c, d) Same for Irene Skylake (48 Skylake processors per node).

Experiments were run on Joliot-Curie, at Très Grand Centre de calcul du CEA (the French Alternative Energies and Atomic Energy Commission (CEA) – TGCC).

Considering that FALL3D-8.1 uses the Message Passing Interface (MPI) for 3D domain decomposition with freedom for the user to choose the number of processors along each spatial direction, to identify the optimal running configuration on Irene ROME and Irene Skylake, a few benchmark cases (with a grid size similar to that of the real benchmark ones, ∼50 million grid points, and 12 particle bins) were ran changing the configuration of nodes and cores used. Results are shown in Fig. C1.

As observed, for this particular grid size, parallel efficiencies are substantially better at Irene ROME, at >90 % with up to 2048 processors (16 nodes). At the Irene Skylake partition, parallel efficiencies already drop to ∼70 % with only 1036 processors (32 nodes). Scalability breaking at a larger number of processors occurs because the number of grid points per subdomain becomes less than the specified range in which communications start to overtake computations (a larger grid size would be needed to sustain speedup ratios that are close to optimal above 2048 processors). Then, considering the resolution of our domain (0.025^∘) and the total grid points of ∼35 million ($\mathrm{1040}\times \mathrm{920}\times \mathrm{35}\times \mathrm{14}$), we fixed the number of nodes to 16 and the number of cores to 768. This configuration allows decomposing the grid points into $\mathrm{32}\times \mathrm{24}\times \mathrm{1}$ (X, Y, Z) subdomains of more than 30 points per spatial dimension. As a result, we increased the speedup 16-fold and the parallel efficiency was fixed to ∼90 %.

Appendix D: Complementary maps

D1 Exceedance probability maps at FL250

Figure D1 shows the exceedance probability maps computed at FL250.

Figure D1 Exceedance probability maps at FL250: the isolines show the probability of reaching an ash concentration above 0.2 mg m⁻³ (a, d), 2 mg m⁻³ (b, e), and 4 mg m⁻³ (c, f) at FL250 at some time from the onset of the eruption to up to 48 h after its end.

D2 Persistence maps at FL250

Figures D2 and D3 show the exceedance probability maps computed at FL250.

Figure D2 Persistence maps at FL250 (large class): the isolines show the probability of reaching or exceeding an ash concentration above 2 mg m⁻³ at FL250 at 1, 3, 6, 12, 18, and 24 h during the eruption to up to 48 h after its end.

Figure D3 Persistence maps at FL250 (medium class): the isolines show the probability of reaching or exceeding an ash concentration above 2 mg m⁻³ at FL250 at 1, 3, 6, 12, 18, and 24 h during the eruption to up to 48 h after its end.