Assessing potential impact of explosive volcanic eruptions from Jan Mayen Island (Norway) on aviation in the North Atlantic

Volcanic eruptions are amongst the most jeopardizing natural events due to their potential impacts on life, assets, and environment. In particular, atmospheric dispersal of volcanic tephra and aerosols during the 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 assess5 ments on the possible consequences for air traffic of a future ash-forming eruption are nonexistent. This study presents the first comprehensive long-term volcanic hazard assessment for Jan Mayen volcanic island in terms of ash dispersal and airborne tephra concentration at different flight levels. In order to delve in the characterization and modelling of that potential impact, a probabilistic approach based on merging a large number of numerical simulations is adopted, varying the volcano’s Eruptive Source Parameters (ESPs) and meteorological scenario. Each ESP value is randomly sampled following a continuous Proba10 bility Density Function (PDF) defined from the Jan Mayen geological record. Over 20 years of climatic data are considered in order to explore the natural variability associated with meteorological conditions and 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 at different flight levels in the atmosphere. The resulting maps represent an aid to civil protection, decision makers and aviation stakeholders in assessing and preventing the potential impact 15 from a future eruption at Jan Mayen.

. Overview map (a) of the study area with the location of structural elements identified on potential field data. Structural elements map (b) for the JMMC: mapped faults, fractures zones and lineaments based on (Peron-Pinvidic et al., 2012;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). Table 1. Possible relative eruption scenarios on Jan Mayen Island. The categorization is based on the volume of tephra emitted in DRE. Data obtained from http://icelandicvolcanoes.is/?volcano=BEE, (Gjerløw et al., 2016). According to the geological record (extending beyond Holocene), subPlinian/Plinian events are highy unlikely (1%). Because of this, they are not included in this is also difficult to assess due to erosion and superimposition of newer vents (possibly covering and removing older ones).
However, considering visible evidence, the (under) estimated number eruptions over the last 10k years is around 45, resulting 90 in an eruption frequency of 1 eruption every 220 years. The duration of the eruptions from Sør-and Midt-Jan is still unknown.
The unrest episode recorded in 1732 (Eggoya, Midt-Jan), which led to the largest known explosive eruption, was a Surtseyan eruption that dispersed tephra over large parts of Jan Mayen and the surrounding seas. The volume of tephra ranges between 0.3-0.4 km 3 .
-Small eruptions are mostly effusive events characterised by small lava flows or small scoria cones, with erupted volumes 100 ranging 10 7 -10 8 m 3 (total less than 0.1 km 3 DRE), corresponding to eruption magnitudes 1 to 2, hence VEI=2. Based on historical occurrence, this scenario can last for about 35-40 hours.
-Moderate eruptions include subaerial, sub-glacial and even surtseyan eruptions depending on 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, moderate eruptions are characterised by aa-lavas but also pahoehoe-105 flows. Surtseyan eruptions are expected to be located on Jan Mayen and 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 days, generate a volcanic plumes between 3 and 11 km in high, and have a range of total erupted volume of 10 8 -10 8.7 m 3 (total volume emitted would be between 0.1 and 0.5 km 3 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 days and one 110 month. As a result, ash fall phases are expected, producing deposits more than one-meter thick within 5 km from the vent. The reference eruption for the Surtseyan type is the Eggoya 1732 AD eruption that produced at least 0.3-0.4 km 3 of tephra (0.16-0.21 km 3 DRE) (Gjerløw et al., 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 characterised by plume heights reaching more than 115 10 km and a range of total erupted volume of 10 8.7˘1 0 9 m 3 (total volume emitted > 0.5 km 3 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 3 DRE (Siggerud, 1972). As the eruption proceeds it becomes more effusive lasting for 1-4 days.

120
-Sub-Plinian eruptions include sub-Plinian to Plinian eruptions characterised by column heights from 15 km to 25 km and a range of total erupted volume of 10 9˘1 0 9.7 m 3 , hence eruption magnitudes of 5 to 6, corresponding to sub-Plinian type I or VEI ≥5. According to (Gjerløw et al., 2016), in the geological record (extending beyond Holocene) there is evidence of 10 tephra layers from subPlinian/Plinian events in 119k years. Because of this, we assign a subjective probability of 1% to this category in case of eruption.

Probabilistic hazard assessment approach
Until a few years ago, volcanic hazard assessment was largely based on the concept of "eruptive scenario", characterized by subjectively-defined eruption conditions. Hazard was then quantified under the strong assumption that the next eruption from a given volcano will be similar to the selected "representative eruptive scenario". However, when assuming a representative eruptive scenario, one is implicitly neglecting the large uncertainties (both aleatory and epistemic) in the parameters that define 130 the scenario also called "intra-size variability".
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" (e.g. Bonadonna et al., 2005). However, the use of a specific and limited range of eruption parameters continues introducing a large biased and uncertainties in the description of potential eruptive processes. For this reason, more recent approaches are based 135 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 (Sandri et al., 2016). The natural approach to do it is by defining probability distributions within plausible ranges. As a consequence, eruption parameters (e.g, total erupted mass, duration of the fallout phase, mass eruption rate, total grain size distribution, etc.) are defined and randomly sampled from specific probability dis-tributions (Sandri et al., 2016). The processes for sampling and weighting possible statistical combinations of values for the 140 volcanological parameters corresponds 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): The four colors cover the erupted volume ranges in the four "classical" eruption categories for Jan Mayen, usually synthesized in 4 representative scenarios with a fixed mass, neglecting the variability in volume around these scenarios. The area under the different parts of the plot correspond to the probability of an Effusive, Medium, High and Subplinian category range eruption respectively, conditional to eruption occurrence. These values are in agreement with previous studies for Jan Mayen (Larsen et al., 2017;Gjerløw et al., 2016).
1. A very broad range of possible eruptive categories, characterized by the total erupted volume, is selected as explained in Section 3.1.
2. The total erupted volume is used to define the total erupted mass, the eruption magnitude, and the VEI. 3. The eruptive category range is split into eruption classes linked to representative members (see section 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 7 − 10 10 m 3 ), up to 6 different truncated Probability Density Functions (PDF) are tested to describe the conditional probability of these 4 mutually exclusive categories: Normal, Exponential, LogLogistic, LogNormal, Gamma and Weibull. The best model is selected according to the Akaike 150 Information Criteria. Indeed, the assumption of a common PDF for the total erupted volume across the different eruption classes allows a smooth and coherent linking among them (Sandri et al., 2016). For JM, the Weibull PDF (see Figure 2) better fitted the expected frequencies on the sub-ranges for the 4 different eruption classes. This PDF is used to assign a conditional probability of occurrence to each simulation as a function of the associated total erupted volume.
5. Considering the behavior of similar scenarios including wet plumes, for Medium and Large classes we account for 155 particle aggregates assuming different aggregate bin. They are characterised by densities in the range of 250 and 350 kg/m 3 and diameters between 100 and 250 µm.

Pulsating eruptions: Modelling and strategy
A novel strategy is proposed to describe the styles and model dispersal from pulsating eruptions, characterized by a series of discrete short-lived events followed by occasional interruption of the emission of tephra. The strategy has been developed 160 considering the ranges of all the ESP described in section 3.1. For each pulsating scenario, the ESP associated with 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 Figure 3): Figure 3. 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.
1. Sampling randomly both the total erupted mass (TEM T ) and total duration of the eruption (Dur T ) considering values reflected in Section 3.1.

165
2. If the sum of masses erupted by all pulses does not equal or exceed the total erupted mass previously sampled, loop to: -Create the i − th 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 lower limit 3 km, peak at 6 km, and upper limit 11 km.
-Compute the total erupted mass for such pulse (TEM i ) using the (Mastin et al., 2009) relationship. (TEM i )), being n the number of pulses generated so far: thereby avoiding small pulses. Compute the new column height (H i ) using (Mastin et al., 2009).
-Else, if total erupted mass obtained ( n i=1 (TEM i )) < 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 ( generate n-1 inter-pulses at random time (Res i ) so that their sum equals δ (δ=Dur T − n i=1 (Dur i )) and insert them (Dur i ). This case actually supposes a continuous eruption where each pulse occurs without a rest period.

Vent location sensitivity
Given the scales of JM Island and the considered domain, the effects of the uncertain vent location on the resulting long-range 180 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 2 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 Figure 4 inset). At these locations we inspected: -Local wind profiles: Figure 5 shows vertical profiles of wind speed and direction averaged for the whole month of December, 2019. As observed, there are little differences in patterns between the two locations.

185
-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.
As expected, the location of potential JM vents does not influence on the ash dispersal pattern. As a result, we will not consider the uncertainty on the vent location and assume a fixed vent at the middle of the island.

Hazard maps and uncertainty quantification
Hazard maps and probability maps (Elefante et al., 2010) are powerful tools to provide information on 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 195 phenomenon (e.g., tephra load at ground level or airborne tephra concentration at different flight levels).
In that regard, our objective is to show the usefulness of HPC-PVHA (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 UK, see Figure 4) from a potential eruption at Jan Mayen Island, while quantifying how the ESP and wind patterns (velocity and direction) influence hazard and probability maps of ash dispersal and airborne tephra 200 concentration.
In this study we try to answer the following questions: which is the probability that, in case of an eruption in Jan Mayen, the ash cloud concentration will exceed the critical condition for safe flights within a domain extending down to the UK airspace after 3, 6, 12, 24 hours since the beginning of the eruption?

205
in case of an eruption in Jan Mayen, which is the probability that airports in Iceland and UK will be affected by the presence of ash?
which is the probability to exceed a predefined hazardous temporal persistence of unsafe flight conditions?
which flight level is likely to be predominantly affected by critical concentrations of volcanic ash?
To this end, considering the impacts of volcanic ash on jet engines summarized in Figure 7, we analyze the results using isolines   -Exceedance probability maps: reporting the probability of reaching ash concentration above a given threshold (0.2, 2, 215 and 4 mg/m 3 ) at different flight levels and anytime during the eruption up to 48 hours 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 3 ) with a 5% probability. has been subjectively selected. However we highlight that our method allows a potential end user to explore any value of exceedance probability: here, for the sake of brevity, we only show the 5% maps as an example.  probability of 5% between 0 and 48 hours after eruption. Black circles correspond to the airports previously cited. (Kristiansen et al., 2012) have 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   Figure 13 shows, from top-left to bottom-right, the probability of reaching or exceeding ash concentration above 2 mg/m 3 at 5000 feet elevation for more than 1, 3, 6, 12, 18 and 24 hours, respectively, from the onset of a large eruption up to 48 hours 245 after its end. Figure 14 shows the same but for the medium-size-eruption and, in the appendix, Figures D3 and D4 display the same information as Figures 13 and 14 but for 25000 feet.

Discussion
In this section, we will discuss the results obtained for each type of analysis proposed in section 4.1.  probability of 5% between 0 and 48 hours after beginning of the eruption, for large and medium eruptions respectively. Among the most important conclusions, we can highlight that an ash-rich eruption originating from Jan Mayen volcano has potential to affect the air traffic over Iceland (after 36 hours) and, to some extent, the Faroe Island, after 48 hours.
255 Figure 15 shows the evolution of the probability over time of having exceeded ash concentration of 2 mg/m 3 at 5000 feet and at international airports of Keflavik and Akureyri (Iceland), Vágar (Faroe Islands) and Luton and Heathrow (UK, London) since the beginning of the eruption. We can observe that the probability of exceeding the threshold at any airport is almost zero during the first hours (ten hours for the medium size and 15 hours for the large case, approximately) and increases until    These results are in agreement with the results shown in Figures 11 and 12, where different maps, at different levels of confidence, produced by cutting the hazard curves at different percentiles depict the relative uncertainties related with airborne ash cloud concentrations and extent for both eruption categories at 5000 feet.
When analyzing the extent, the influence of the eruption category is also relatively important. At 5000 feet, the volcanic ash 280 cloud with concentrations up to 2 mg/m 3 , affecting maximum risk actions like takeoff or landing, could reach a large part of the UK with a probability between 2.5 and 50% for both eruption classes. It could threaten the vast majority of flights to and from northern routes. However, at 25000 feet, the ash cloud would affect almost the entire FIR, even reaching part of the UK for concentrations up to 0.2 mg/m 3 when large size eruptions. For medium-size eruptive class, only polar routes above 25000 feet would be threatened. Then, we can conclude that for medium-size eruption class, only polar routes above 25000 feet would be 285 threatened.
Finally, similar to other types of analysis such as ground load of tephra and Probabilistic Seismic Hazard Assessment, Figure   16 provides a graphical representation of relative uncertainties related with airborne ash cloud concentrations above 0.2 mg/m 3 , 2 mg/m 3 , and 4 mg/m 3 at 5000 feet at Keflavik airport. This result, that we quantify at each point of the target domain, allows integrating hazard in quantitative risk analysis, through fragility curves. In this view, it represent the most complete way to 290 quantify hazard. Specifically, no dramatic differences are seen depending on the eruption size, and there is a non-negligible probability to overcome the 2 mg/m 3 threshold, even for low percentiles, given an eruption.

Persistence maps
According to Figure 7, engines exposed to such persisting concentration conditions for more than 3 hours would require 295 inspection. Therefore, we conclude that in terms of ash extent, results (Figures 13 and 14) are slightly different: concerning altitude, at 5000 feet a large part of FIR (reaching in some extent UK) with probabilities 5 and 50% would be affected for both eruption classes, while at 25000 feet ( Figures D3 and D4)), such conditions would affect only high latitude air routes (above 68°N).
However, when analyzing the spatial pattern for long-term persistence (more than 12 hours), we find some differences 300 depending on the eruption category. For persistence above 12 hours, at 5000 feet, an ash cloud from a large/medium eruption has 2% to 10% probability to reach latitudes as low as 65°N/ 62°N. Such southernmost latitude increases for longer persistence values, meaning that (obviously) only closer to the source we may get long-persisting clouds.
Finally, Figure 17 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 at different airports in Iceland and UK for at least 24 hours since the beginning of the eruption. The most affected airports are Akureyri, Vagar and Keflavik (London LTN and London LHR have very low probabilities (1.5%-2%) of persistence of 1 to 6 hours only when the simulation is carried out for 48 hours after the beginning of the eruption). Concerning the level of persistence, we can highlight that both eruptions categories have similar behavior. Scenarios with persistence greater than 18 hours are highly unlikely.
However, when analyzing probabilities, medium eruptions reach persistence probabilities twice higher than large ones. This 310 observation can again be associated with their eruptive dynamic. The sustained injection of tephra into the atmosphere related with a series of discrete short-lived events increases the probability of prolonged persistence scenarios. Although limitation due to the lack of a complete geological record composed of both chronological and statistical data can be a bias for hazard assessment on Jan Mayen Island, the spatio-temporal analysis and threat assessment presented here do provide 315 the first comprehensive analysis of potential impact on aviation safety in the north Atlantic due to a future explosive activity.
We believe that, despite the intrinsic limitations of the methodology (partially due to the scarcity of data), this work represents an important element in the long-term volcanic hazard assessment for Jan Mayen volcano of interest for key stakeholders like aviation authorities, VAACs and governments.
Among the most important conclusions, it can be highlighted that, an ash-rich eruption originating from Jan Mayen volcano 320 has potential to affect the air traffic over Iceland (after 36 hours) and, to some extent, the Faroe Island, after 48 hours. Concerning airborne ash concentration and extent, for large eruptions, concentrations above 2 mg/m 3 (even 4 mg/m 3 , originally considered no fly zone) would affect part of the Icelandic airspace (at different flight levels) with exceedance probabilities between 5 and 50% at some time during the eruption up to 48 hours after its end. For medium eruptions, these dangerous concentrations would affect only low flight levels. Above 25000 feet only polar routes would be affected.

325
When analyzing persisting concentration conditions where aircraft engines are exposed to high concentration for more than 3 hours, we conclude that at 5000 feet a large part of FIR (reaching in some extent UK) with probabilities 5 and 50% would be affected for both medium and large eruption categories. At 25000 feet 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, that should be routinely 330 considered in all this kind of studies.
Code and data availability. Scripts and pipeline programs will be uploaded according to the manuscript preparation guidelines This section provides an overview about the High Performance Computing (HPC) environment used in this study and the setup process associated with the Fall3d model to simulate the eruptive scenarios. We will describe the most relevant settings 335 to optimize both energy and computational resources, as well as the simulation scheme followed.

Appendix A: Simulation setup
In order to meet computational constraints, we run the experiments by months. We randomly sampled 1500 wind fields on the time interval 1999-2019 from ECMWF Reanalysis database for each eruptive class (large and medium). Then, we run 1500 simulations (scenarios) combining meteorological conditions and volcanological parameters for each category. Since medium 340 eruptive classes are characterized by a series of discrete short-lived events, the total number of scenarios for such category was 3763.
As can be seen in Figure A1, which shows the variance of the tephra concentration for a given grid-point with respect to the number of scenarios simulated, from 900 scenarios, the variance of the concentration begins to stabilize. This stabilization also suggests a reduction in uncertainty related to the intra-size variability of the eruptive scenarios themselves. To run the 345 scenarios, we used the Fall3D-8.1 model (Folch et al., 2020), an open-source off-line Eulerian model for atmospheric passive transport and deposition which solves a set of advection-diffusion-sedimentation (ADS) equations on a structured terrainfollowing grid using a second-order Finite Differences (FD) explicit scheme. The simulation scheme was run on a (approx.) Large eruptive categories. Figure A1. Value of the variance of tephra concentration for a given grid-point with respect to the number of scenarios simulated.
Appendix B: Sampling and processing workflow As mentioned above, in order to achieve an unbiased PVHA, we should account for the uncertainty in the ESP, as well as the magnitude and intensity variability in the eruption class, in the vent position, and in the wind field. The natural approach to do it is by simulating and appropriately combining a statistical number of representative eruptive scenarios, where each 355 of these representative eruptive scenarios is characterized by randomly-defined eruption conditions. Such eruption conditions are composed by several eruptive parameters. Then, for each eruptive class, we set the PDF for each eruptive parameter as in (Sandri et al., 2016). It is important to note that in this study we only addressed the medium and large eruptive classes. Table   B1 summarizes the PDFs and values ranges of the main eruptive parameters for Jan Mayen island. distribution and sphericity of tephra particles from their PDFs.
2. Compute the mass fraction (%) associated to tephra fallout with respect to the total erupted mass according to the available estimations from field data analysis. For medium eruptive class, the value of % is fixed to 0.8, meanwhile for large eruptive class is randomly sampled from [0.05, 0.10].
4. Sample a time for the eruption start over a period of 20 years (1999-2019) considering the corresponding meteorological fields for the duration of the fallout phase, and associate this randomly to a combination of the volcanological parameters.
For this, we download ECMWF ERA5 reanalysis meteorological data associated with the date and duration of the eruption.

370
5. Run FALL3D to obtain the tephra loading at different flight levels. Modify FALL3D's input file with both, meteo and newly sampled data values. If the eruptive scenario assumes tephra composed of more than one type of particle, it will be necessary to create as many input files as types of particles and modify the source block independently for each of them 6. Once the eruptive scenario has been built, it is run in the framework of High performance Computing (HPC). To do that, 375 we use TGCC Joliot Curie, the most powerful French supercomputer dedicated to French and European research.
7. The probability of each combination is weighted in accord with the associated magnitude. The results obtained are computed using PVHA-WF software. As a result, we obtain hazard and probabilistic maps describing the airborne ash concentration and time-persistence at different flight levels on a large-scale and high-resolution domain.
Typical tephra particle densities and total grain size distributions were chosen consistent with previous values reported for 380 Eggoya (1790) and Grimsvotn (2004) eruptions for Medium and Large categories respectively.
Regarding the type and spatiotemporal distribution of aggregates (ash particles with diameters <63 µm settle from eruption clouds and adhere to each other as a result of electrostatic attraction, moist adhesion between particles and hydrometeor formation) were also chosen consistent with previous values reported for similar surtseyan and phreatomagmatic eruptions.

385
Experiments were run on Joliot-Curie, at CEA/TGCC (a relevant HPC environment dedicated to French and European research). We had two partitions assigned, Irene-rome and Irene-skylake (see Table C1 for specifications) with a total of 63M core hours.
When working in HPC environments, it is mandatory to find the optimal configuration of both the number of cores and Table C1. Joliot-Curie supercomputer. Characteristics corresponding with the two partitions available on this study.
Irene-Skylake CEA/TGCC 1656 Intel-skylake 2.7GHz bi-processor compute nodes with 48 cores per node (24x2).This totals 79488 cores with 180GBy/node. nodes used in order to optimize the energy consumption and computing time. Considering that FALL3D-8.1 uses MPI for 3D 390 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, we run a few benchmark cases (with grid size similar to that of the real benchmark ones, 50M grid points and 12 particle bins) changing the configuration of nodes and cores used.
Results are shown in Figure C1.
As observed, for this particular grid size, parallel efficiencies are substantially better at Irene-rome, with >90% up to 2048 395 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 sub-domain becomes less than the specific range in which communications start to overtake computations (a larger grid size would be needed to sustain speed-up ratios close to optimal above 2048 processors). Then, considering the resolution of our domain (0.025º), and the total grid points 35M (1040*920*35), we fixed the number of nodes to 16 and the number of cores to 768.

400
This configuration allows decomposing the grid points into 32*24*1 (X,Y ,Z) subdomains of more than 30 points per spatial dimension. As a result, we increased the speed-up 16 times and the parallel efficiency was fixed to 90%. Mingari contributed to coding the scripts and software, analysis and/or interpretation of data and drafting the manuscript.