Modelling pyroclastic density currents of the April 2021 La Soufrière, St Vincent eruption: from rapid invasion maps to field-constrained numerical simulations Open Access

After 41 years of quiescence, 2 months of minimal seismic unrest and the extrusion of a lava dome in the crater that started in late December 2020, La Soufrière, St Vincent erupted explosively on 9 April 2021, and entered a new eruptive phase that lasted 2 weeks. The response team, assembled at the time of the first extrusive activity, successfully forecast the transition to explosive activity, enabling evacuation of areas at the highest risk around La Soufrière on 8 April. Thanks to this successful crisis management, no fatalities occurred during this eruption. Part of this success was due to the clear identification of the potential volcanic phenomena that could result from the eruption of La Soufrière and the destruction of the newly formed lava dome. Among them, pyroclastic density currents (PDCs) were considered a serious threat. These currents comprise a mixture of hot volcanic particles and gases that can form by a variety of mechanisms including lava dome collapse, lava dome explosion or eruption column collapse. Destructive and fatal PDCs had already occurred during the 1902–03 (Anderson and Flett 1903) and 1979 eruptions (Shepherd et al. 1979) of La Soufrière, as well as from volcanoes on other islands in the Caribbean such as Mount Pelée in Martinique in 1902–03 (Lacroix 1904), or Soufrière Hills Volcano in Montserrat from 1995 to 2010 (Druitt and Kokelaar 2002; Wadge et al. 2014). From the Integrated Hazards map of St Vincent (modified from Robertson 2005; see also Fig. 1b), PDC hazards are considered together with all other volcanic hazards in the ‘very high hazard’ zone that corresponds roughly to the extent of the La Soufrière volcano edifice. However, two towns are close to the edge of this zone: Georgetown, on the eastern coast, lies just inside the boundary of the high hazard zone and Chateaubelair, along the western coast, lies outside the high hazard zone. A serious concern therefore is the possibility of unpredicted, voluminous PDCs capable of inundating a particularly large and densely populated area. In order to correctly assess the PDC hazards during the 2021 eruption, the response team decided to apply numerical modelling of PDCs to try to answer two main questions raised at the beginning of the April crisis: ‘Which areas of La Soufrière volcano are the most likely to be affected by PDCs?’, and ‘What parameters would generate a PDC that inundates Georgetown or Chateaubelair?’. A new PDC modelling project was then quickly set up to try to answer those questions during the crisis.

A typical way of assessing complex PDC hazards is to use numerical modelling (Calder et al. 2015). The outputs can then be integrated with other volcanic hazards such as lahars or tephra fallout in volcanic hazard maps (Neri et al. 2015; Newhall and Pallister 2015). There are currently two main approaches to generating volcanic hazard maps (Calder et al. 2015; Rouwet et al. 2019): (1) a deterministic, or scenario-based approach, largely used in the past 40–50 years, based on several scenarios representative of the volcanic historical records or used as proxies for possible future activity; and (2) a probabilistic approach, increasingly used in the last decade, that relies on the statistics of recurrence or an aleatory exploration of input parameters (Tadini et al. 2017; Sandri et al. 2018; Tierz et al. 2018). The probabilistic approach does not require the calibration of input parameters prior to simulations, as a range of values is explored during the modelling process instead, but it requires a significant time to complete all simulations, and considerable computational capabilities. Thus, the probabilistic approach is well adapted for pre-eruptive and long-term hazard assessment. The deterministic approach on the other hand allows for a faster forecast as fewer simulations are needed, making it more suitable for short- to medium-term forecasting. A drawback, however, is that it does require the user to define input parameter values, for which uncertainties are usually difficult to estimate (Rouwet et al. 2019). Recent advances have, however, shown significant improvements in the estimation of input and output model uncertainties using Bayesian metrics (Charbonnier et al. 2018; Gueugneau et al. 2020; Mead et al. 2021). Because of their respective advantages and drawbacks, combining deterministic and probabilistic approaches is now considered more beneficial than choosing only one of them, as it guarantees a more complete hazard forecasting at many levels (Newhall and Pallister 2015; Rouwet et al. 2019). Regarding La Soufrière, St Vincent, probabilistic hazard assessment has only been recently investigated by the University of the West Indies Seismic Research Centre (UWI-SRC) (when the first unrest appeared), and the use of a complementary deterministic approach to investigate the two main questions raised above appeared to be timely and suitable in the case of the 2021 eruption. It allowed us to test specific scenarios in a short timescale, adapted for such a syn-crisis hazard assessment.

In this study, we present results of numerical simulations aimed at (1) assessing PDC hazards during and after the April 2021 eruption, and (2) better constraining input parameters of numerical simulations to improve rapid syn-crisis hazard assessments generally. The study is therefore organized in two steps: (1) a syn-crisis phase, in April 2021, during which PDC hazards were assessed with the aim of supporting the UWI-SRC by making rapid invasion maps; and (2) a post-crisis phase during which syn-crisis maps were reanalysed in the light of new data from fieldwork. For the first phase, we ran numerical simulations of PDCs at La Soufrière using the two-phase version of VolcFlow, during the 2 weeks after the onset of explosive activity on 9 April. We chose this code for: (1) its capacity to simulate the extent of both dilute and concentrated regimes of PDCs during a single simulation; (2) its rapidity, i.e. simulation run <5 h; and (3) the simplicity of generating invasion maps of concentrated and/or dilute PDCs from simulation outputs. Cumulative PDC invasion maps were created for two sectors of La Soufrière, i.e. the southeastern and southwestern flanks. In the second phase, thanks to new post-eruption data as a PDC deposit map, new simulations were run and the best-fit simulation with respect to the deposit map is presented. Cumulative invasion maps and the best-fit simulation are compared, and the differences are discussed. Misfits in values of input parameters, differences in data processing, the capability of the two-phase version of VolcFlow for rapid hazard assessment, and the pros and cons of a scenario-based approach for PDC hazard assessment are discussed in terms of lessons learned for future eruptions.

La Soufrière volcano is a 1220 m high stratovolcano located in northern St Vincent Island, in the south of the Lesser Antilles volcanic arc in the Caribbean Sea (Fig. 1a). Several eruptions have been recorded in historical times (1718, 1812, 1902–03, 1971–72, 1979) making this one of the most active volcanoes in the Lesser Antilles chain. These eruptions, which have comprised both effusive and explosive phases, have impacted both the populations of St Vincent and the wider region. Scoria-rich PDCs have been mapped from several of these eruptions (1718, 1812, 1902–03, 1979; Cole et al. 2019). The major eruption of 1902–03 (VEI 5, Cole et al. 2019), which generated extensive PDCs, resulted in c. 1500 fatalities and a considerable economic impact. The last eruption in 1979 produced a series of 13 explosions and resumed with a lava dome growth at the centre (Shepherd et al. 1979). Since then, the volcanic activity reduced to a fumarolic field along the southern edge of the 1979 lava dome.

In November 2020 an increase in background seismicity at La Soufrière was observed and site visits indicated minor changes to fumarolic activity on the 1979 dome and the small crater lake (Joseph et al. 2022). On 27 December 2020 the NASA Fire Information for Resource Management System detected a thermal anomaly inside the summit crater, which was confirmed as new dome growth by visual observations on 29 December 2020. This effusive activity continued for a period of 3 months with a steady rate of dome growth until rapid extrusion, increased venting and visible incandescence were observed (6–8 April 2021), accompanied by the appearance of a banded tremor of increasing magnitude on 8 April 2021 (Joseph et al. 2022). At this time the alert level was raised to red, triggering the evacuation of approximately 16 000 persons from the highest risk areas. Explosive activity initiated on 9 April at 12:41 UTC and consisted of 32 discrete events with plumes up to 15+ km high (Joseph et al. 2022). As the activity intensified, the visual identification of discrete events became more difficult but was confirmed using seismicity and satellite observations. Later explosive activity of mainly vulcanian style lasted until 22 April 2021 (Joseph et al. 2022; Miller et al. 2022), forming a new crater inside the pre-existing La Soufrière crater (Fig. 1c), modifying the morphology of the summit area with the enlargement of the new crater, and an important accumulation of fresh pyroclastic material around it within the pre-existing crater (Fig. 2b). Following this explosive phase was a period of gradually declining activity, with daily seismicity and sulfur dioxide flux returning to close to background levels in March 2022.

The PDCs were associated with the explosive activity and emplaced mainly on the southwestern flanks of the volcano. The PDC deposit map created after field campaigns in early 2022 is presented in Figure 1c (see also Cole et al. 2023). For context, this section briefly summarizes PDC characteristics, and the reader should refer to Cole et al. (2023) for a more detailed description of the PDC deposits and associated eruptive stratigraphy. Note that data presented in this section were gathered mostly after the eruption, and we only had access to the position of the crater (Fig. 2a) during the syn-crisis phase.

A precise timeline of all PDC events could not be created owing to the lack of visual observations. For the first 48 h of the eruption, copious ash in the atmosphere as a result of closely spaced explosions caused darkness, making visual observations difficult or impossible. Satellite imagery (radar backscattering, Fig. 2c) indicated that the first PDCs had been emplaced by 11 April 2021 in the southern and western sectors, extending 2–2.5 km from the crater rim. Based on stratigraphic studies and eyewitness accounts, Cole et al. (2023) infer that the first PDCs probably occurred late on 10 April. The first visual observations of PDCs were made on 12 April 2021, at 08:15 UTC, with flows extending down multiple valleys in the southern and western sectors of the volcano. PDC deposits were visually confirmed, via boat, having reached the sea in the Larikai and Roseau valleys on 12 April (Fig. 2c). Fieldwork and satellite imagery showed that some minor emplacement of PDC deposits occurred in the southeastern sector of the volcano close to the crater rim and extended short distances (c. 1 km) down valleys to the SE (effectively at the head of the Rabacca valley; 14 April 2021 at 02:48 UTC). They also extended down several of the channels leading from the southern part of the crater rim and into the upper Wallibou valley, reaching about as far as Trinity falls/Wallibou hot springs, approximately 2.5 km from the crater rim.

Deposits of the valley-filling PDCs were studied in both the Roseau and Larikai valleys (dark red areas, Fig. 1c). They display a succession of three to five PDC units of 1–4 m thick, giving a total deposit thickness of 20–30 m.

Dilute PDCs, or ash-cloud surge, were extensive in proximal regions (Fig. 1c, light red area), but had slightly shorter runouts than the valley-filling PDCs, extending generally 1.5–2 km from the crater rim. These dilute PDCs caused extensive tree damage in the inundated area, above the valley thalwegs (on slopes and ridges), and emplaced thin deposits, <15 cm thick, made from fine lapilli, typically poor in fine ash. These deposits are relatively continuous but do show some minor lateral thickness variations.

We used the two-phase version of the numerical code VolcFlow developed to simulate the extent of both regimes in pyroclastic currents, i.e. concentrated and dilute (Kelfoun 2017). The code is based on two coupled, depth-averaged currents: one for the basal concentrated zone and one for the overriding dilute zone – the ash-cloud surge. The dynamics of each current are modelled using depth-averaged equations of mass and momentum balance in the x and y directions. The ash-cloud surge requires an additional equation, as its density varies in time and space owing to the loss of mass through sedimentation. The two layers are then coupled and exchange mass and momentum following two exchanges laws. To simulate stresses applied to the concentrated layer during transport using a depth-averaged approach, the plastic rheological law is used, involving a constant retarding stress T, which has been shown to successfully reproduce various features of such currents and their deposits (Kelfoun et al. 2009, 2017; Kelfoun 2011; Charbonnier and Gertisser 2012; Ogburn and Calder 2017; Gueugneau et al. 2019, 2020). The ash-cloud surge is simulated as a turbulent continuum that loses momentum owing to turbulent drag stresses. The complete description of the physical model, the equations and all the parameters used in VolcFlow are summarized in Kelfoun (2017) and Gueugneau et al. (2019; 2020, supplementary materials).

Simulations were undertaken using the 2 m resolution Pleiades digital elevation model (DEM) (freely available on 7 April – R. Grandin, IPGP), downgraded to 8 m resolution to allow faster simulation runs. Owing to limitations on the southern extent of this DEM, we merged this DEM with a 15 m resolution DEM covering the entire island, published by UWI-SRC in 2005 (Robertson 2005). The merged DEM was then resized to 8 m cells for consistency with the first simulations.

We ran preliminary simulations on 10 April to quickly identify sectors of the volcano most at risk of inundation. Based on these results, the simulation domain could then be reduced to save computation time for the subsequent simulations. Thanks to the radar images acquired on 10 April (Fig. 2a), the position of the new crater was located, and test simulations could be run with a source centred in this new crater. As the exact source mechanism for the PDCs was not identified at this time, e.g. either a column or a lava dome collapse, our simulated source conditions were simplified to a constant mass flux supplied during 75 s from a circular shape with a diameter that encompasses the entire crater, as well as the external side of the crater rim (without any topographic correction; see Fig. 3). We note that such source conditions can be applied to simulate both gravitational dome and column collapses and that we disregarded any possible syn-eruptive topographic changes that could have occurred inside the summit area at that time. Owing to prior experience acquired using VolcFlow on three different PDC events elsewhere (Kelfoun et al. 2017; Gueugneau et al. 2019, 2020), the values of the remaining input parameters were estimated and calibrated accordingly (Table 1). Using such source conditions, the total PDC volume was purposely overestimated (50 Mm³) to encompass all probable scenarios considered at that time. Figure 3 shows the simulated ash-cloud surge (in transparent blue) and concentrated flows (in red) that emerged from the crater and propagated to the south and west towards the sea, channelized in deep valleys (Fig. 3b–e). The simulations showed that PDCs can be emplaced on the southeastern flank of La Soufrière only if they emerge from the south or southeastern lip of the crater. They also showed that only the southern part of La Soufrière could be reached by a PDC owing to major topographic barriers preventing their northwards propagation.

These conclusions and our maps in Figure 3 were shared with UWI-SRC, and two different sectors were selected for further investigation: (1) the southwestern flank, with the highest probability of inundation based on the preliminary simulation; and (2) the southeastern flank. We also noticed that our simplistic source could be refined since most of the material supplied at the source remains inside the crater, as shown in Figure 3. In the rest of this study, this unnecessary supply is avoided by setting a source along the external side of crater rim, with no initial velocity.

To generate an inundation map, a scenario was defined, and a series of simulations was run for which the values of one or several input parameters of the code were varied. To model the complex dynamics of PDCs with the two-layer version of the code VolcFlow, 13 input parameters must be defined prior to each simulation, as detailed in Table 1. Note that the reader is referred to Kelfoun et al. (2017) and Gueugneau et al. (2019, supplementary material) for more information on the influence of each parameter on the model dynamics.

In the case of the April 2021 crisis, none of VolcFlow's input parameters were known prior to the eruption. As for the preliminary simulation, the value of each fixed input parameter was estimated and calibrated from values used in previous studies (Table 1). The last two parameters of the table, i.e. supply duration and total volume, are specific to each eruption and cannot be estimated from previous work. Therefore, we chose to vary the source position/geometry and the total volume to create our scenarios and keep the 12 other parameters fixed in all simulations. Three different volumes were explored to define three PDC-type scenarios:

1. S1, with a volume of 10 Mm³, typical of small-volume PDCs (FlowDat database, Ogburn 2012), similar to the 25 June 1997 PDCs at Soufriere Hills, Montserrat (Loughlin et al. 2002);

2. S2, with a volume of 20 Mm³, twice the targeted volume;

3. S3, with a volume of 30 Mm³, three times the targeted volume, similar to the 8 May 1902 PDCs at Mount Pelée, Martinique (Lacroix 1904; Gueugneau et al. 2020).

Owing to the presence of Georgetown – a significant population centre – and key evacuation routes for settlements on the NE coast, in conjunction with UWI-SRC it was decided to first investigate the southeastern sector on 11–12 April. The southwestern sector was then investigated following this on 13–14 April. The extended zones of Chateaubelair and beyond Georgetown were investigated later, on 16–17 April when we finally acquired the DEM of the whole island to create the merged DEM.

To evaluate the PDC invasion potential on the southeastern sector of La Soufrière, we set an arcuate-shape source along the southeastern quarter of the crater rim (see Fig. 4, black area on cumulative maps). Irrespective of the scenario, using this source, PDCs are channelized into the Waribishy and Rabacca valleys that both drain eastward towards the sea, threatening Orange Hill and minor settlements in this area. Concentrated flows remain confined in these valleys and reach the sea at the mouth of the Rabacca for S2 and S3. An increasing volumetric rate at the source increases the flow runout (Fig. 4), enhancing the capacity of flows to overspill their drainages, especially near the mid-point, where both valleys connect. The simulations indicate that any PDC generated along the southeastern quarter of the crater rim could eventually drain into the Rabacca valley and represent a threat for villages in the vicinity of this river; however, it seems unlikely that concentrated flows would reach Georgetown. Dilute ash-cloud surges covered the entire southeastern flank and spread mostly east towards the sea as Mount Brisbane hills block their southward propagation. Differences in flow extent between the three scenarios are important with an increase in the inundated area by a factor of two between S1 and S3. Although the modelled ash-cloud surge in S3 seems to reach the western edge of Georgetown, the potential impact on the town is inconclusive because the southward extent of the Pléiades DEM is limited. Owing to the limitation of our simulations, including the aforementioned simplified source conditions, the following conclusions were shared with UWI-SRC:

• Areas bordering the Rabacca valley are at risk if a PDC emerges from the southeastern region of the crater rim, as it will eventually flow into that river and potentially reach inhabited areas along the coast.

• Georgetown and inhabited areas along the coast are at risk only in the case of a notably voluminous PDC (e.g. >30 Mm³) that produces an extensive ash-cloud surge. Otherwise, ash-cloud surges seem restricted inshore.

To investigate the southwestern sector, we used the same arcuate-shape source as used previously, but now along the southwestern quarter of the crater rim. As a result, all valleys of the southwestern flank are inundated since they all emerge from the crater rim (Fig. 5). For all scenarios, concentrated flows reach the sea at least at Roseau and Larikai valleys, but also at Wallibou valley for S2 and S3. No concentrated flow is able to escape valleys from that sector and remains well channelized towards the sea. Therefore, Chateaubelair seems out of reach, and only populations around the Wallibou valley mouth are threatened by the passage of a PDC, and only for very voluminous scenarios (S2 and S3; Fig. 5). Here again, the effect of an increasing volumetric rate is mostly restricted to an increase in the flow runout. Dilute ash-cloud surges are extensive in all scenarios and cover the entire southwestern flank, but barely cross Richmond Hill in the south that acts as a barrier. However, similarly to the southeastern sector, the southward extent of the Pléiades DEM prevented us from concluding clearly about Chateaubelair's potential PDC invasion. Owing to the limitation of our simulations, our conclusions shared with UWI-SRC were:

• Concentrated flows remain confined within valleys and do not threaten any inhabited areas, except around the Wallibou valley's sea outlet, for voluminous scenarios only.

• All valleys of the sector can potentially be invaded by a PDC, and the likelihood of a PDC reaching the sea is high, even for small-volume scenarios.

• Chateaubelair seems out of reach of concentrated flow because of the morphology of the topography, but we cannot conclude clearly for the potential inundated area of the ash-cloud surge.

To explore possible ‘worst-case scenarios’, involving potential PDC inundation of Georgetown and Chateubelair, we ran extended sector simulations over the merged DEM (resized into two areas of 8 × 10 and 8 × 8 km for Chateaubelair and Georgetown respectively, Fig. 6). As previous simulations showed that concentrated flows are unable to reach both towns, we only focused on dilute ash-cloud surges propagation. For the Chateaubelair sector (left map in Fig. 6), the source was simplified to a simple spot at the beginning of the upper Wallibou valley, almost perfectly aligned along a straight line from Chateaubelair to the crater, aiming at forcing the surge to follow this direction. However, the topographic control on the surge remains strong and the flow is deflected seaward when it reaches the Wallibou valley. Increasing the volume to 30 Mm³ for S2 yields similar results to those of S3 from the previous study case. We note, however, that the surge begins to get around the Richmond hills towards the coast as the relief decreases and reaches the northern end of Chateaubelair. We then increased the volume further to 50 Mm³ for S3, which caused the surge to extend further south, reaching the coast and inundating two-thirds of Chateaubelair. The radial spreading of the surge flowing over the sea is also enhanced by the artificial morphology of the sea itself, modelled here as a flat surface. This phenomenon has already been visually observed on the south side of La Soufrière of St Vincent in 1902 (Anderson and Flett 1903, see the discussion section), Soufrière Hills at Montserrat in 2003 (Edmonds and Herd 2005), at Krakatoa in Indonesia (Carey et al. 1996) and also numerically at Mount Pelée in Martinique in 1902 (Gueugneau et al. 2020). Although it seemed accurate in Gueugneau et al. (2020), the interaction between an ash-cloud surge and the sea is poorly known and we did not speculate further. The same three scenarios were run on the Georgetown's sector simultaneously. Here the source is an arcuate shape instead of a dot in order to include all tributaries of the Rabacca's drainage system. A smaller arcuate source is used than in Figure 4 to restrict the focus on the Rabacca drainage system only. Similarly, only voluminous scenarios (i.e. >30 Mm³) are able to produce a surge that inundates Georgetown, with an enhanced lateral spreading towards the sea. We note, however, that the Mount Brisbane hills are smoother than the Richmond hills from the SW sector, and oriented in a NE–SW direction rather that east–west, which does not protect Georgetown as effectively as in the case of Chateaubelair. Owing to the limitation of our simulations, the conclusions shared with UWI-SRC were as follows:

• Chateaubelair and Georgetown could be reached by the ash-cloud surge of a notably energetic and voluminous PDC (>30 Mm³), which is not compatible with a Vulcanian or violent Vulcanian-type PDC, but more a blast-like PDC (Mt Pelée in 1902, Montserrat in 1997; Sparks et al. 2002; Gueugneau et al. 2020) or a PDC initiated during a large-magnitude eruption (i.e. Plinian eruption). Therefore, invasion of Georgetown and/or Chateaubelair by a PDC is not impossible, but is unlikely.

In summary, simulations showed that the volcano summit topography has a major influence on determining which sector will be invaded, and that the two sections along the southern half of La Soufrière crater, depicted as arcuate sources in Figures 4 and 5, can be potential emerging points for PDCs. Our maps allowed us to: (1) identify sectors at most risk; (2) identify the worst-case scenarios leading to an invasion of Chateaubelair and Georgetown; and (3) confirm the accuracy of the volcanic hazard map of St Vincent for PDCs, with a warning that Chateaubelair could be inundated by a PDC, but at a low probability. Note that important precautions on the reliability and appropriateness of our results should be considered owing to the uncertainties caused by our lack of knowledge of the PDC properties and ongoing eruptive dynamics. The randomness of the natural system, i.e. the location of the vent, whether a fountain- or column-collapse feeds the PDC, the number of collapse events and PDC pulses, and how voluminous the PDC is, cannot be predicted.

To evaluate our choice of input parameters that were initially selected, using the new data available after the crisis, we ran new simulations using the same DEM and numerical model aiming to reproduce the post-crisis map of PDC deposits as determined by fieldwork and satellite imagery analysis (Cole et al. 2023). We used Bayesian metrics introduced by Charbonnier et al. (2018) and already used for PDCs at Mount Pelée in Gueugneau et al. (2020) to determine the best-fit simulations. These metrics are based on the ratio between the area of the simulated flows (A_sim) and the area of the PDC deposit in the field (A_obs): matching the area between the simulated and observed flows is called ‘true positive’ (TP), the over-simulated area is called ‘false positive’ (FP) and the missing simulated area is called ‘false negative’ (FN). The three coefficients we used were calculated as follows:

• the Jaccard similarity coefficient (R_J) uses TP and the union of areas inundated by observed and simulated flows:

  $RJ=|Aobs∩Asim||Aobs∪Asim|×100=TPTP+FN+FP×100$

  (1)
• the model sensitivity (R_MS) uses TP and the area inundated by simulated flows:

  $RMS=|Aobs∩Asim||Asim|×100=TPTP+FP×100$

  (2)
• the model precision (R_MP) uses TP and the area inundated by observed flows:

  $RMP=|Aobs∩Asim||Aobs|×100=TPTP+FN×100$

  (3)

The result of the best-fit simulation is presented in Figure 7. To help with the comparison, the maximum extent of the real dilute PDC deposit is also reported on the figure as a black line. This shows that the simulation satisfactorily reproduces the real deposit (Fig. 1c), with the same valleys inundated by the simulated concentrated flows, i.e. Larikai, Roseau, upper Rabacca and upper Wallibou valleys, while reaching the sea only in the Larikai and Roseau cases. However, these concentrated flows spread further than the real ones, possibly owing to the low resolution of the DEM (8 m), resulting in relatively low values for the three metrics in the Table 2. In addition, the longer flow runout in upper Wallibou valley did not match the actual PDC. For the simulated ash-cloud surge, its extent is quite close to the observed one in Figure 1c, giving metrics >50% and up to 92.4% for the model sensitivity (Table 2). It does, however, have a shorter runout in the Larikai and upper Rabacca valleys, and a longer runout in the upper Wallibou valley owing to the excess of concentrated PDC volume in this valley.

To match the PDC deposit map, we needed to modify some of our simulation inputs used for the syn-crisis simulations. First, the source was reassessed because no better results were obtained with the single arcuate-shape source used in the syn-crisis simulations. To inundate all valleys of La Soufrière's southern flank, we found that setting four independent arcuate-shape sources at the beginning of the four drainage systems (i.e. Larikai, Roseau, upper Rabacca and upper Wallibou valleys; see yellow shapes in Fig. 7) gave the best results. We recall, however, that even if we used four independent sources around the crater as a better representation of the real source conditions and PDC sequence, the four simulated flows were generated simultaneously at each location (i.e. single PDC pulse). In addition, no bow-shaped source configuration was able to generate an ash-cloud surge that travels far enough eastward and northward. The best explanation for this is that the surge is only partially generated from the concentrated flow, while a part of it is already generated at the source. To test this hypothesis, we added an independent supply of ash-cloud surge, of approximatively the size of the new crater (blue disc in Fig. 7), which acts as a secondary ash-cloud surge source. This source configuration appeared to be the most successful at reproducing the observed surge extent, enabling the dilute current to spread northward and eastward without increasing its runout along the southwestern valleys. Therefore, the surge production coefficient has been drastically decreased by two orders of magnitude to equilibrate the total surge production and reproduce the observed surge extent.

Finally, to accommodate for the smaller runouts and extent of the mapped PDC deposits than those obtained from our simulations, the total volume of the simulation was decreased. The best simulation was obtained for a volume of 5.6 Mm³, half of the volume of S1, the smallest-volume scenario. For finer adjustments, the constant retarding stress of the concentrated flow was increased to obtain shorter runouts, in accordance with previous studies showing that the constant retarding stress is inversely proportional to the flow volume (Kelfoun 2011; Charbonnier et al. 2020).

The ‘post crisis’ best-fit simulation map presented in Figure 7 and in its inputs listed in Table 2 helps to identify weaknesses in the syn-crisis simulation scenarios and the choice of input parameters. The surge production coefficient had the highest variance, with two orders of magnitude difference between syn-crisis and post-crisis simulations and resulting in an overestimation of the ash-cloud surge extent. The value used for the syn-crisis simulation was the highest of the range (see Table 1) used to explore worst case scenarios. In contrast, the best-fit value is in the lowest end of the range, similar to the production coefficients used at Merapi (Kelfoun et al. 2017). This means the values of the surge production coefficient used for this eruption are appropriate and not out of subject, but the spread of the range is an impediment as it can lead to highly variable results. The other noticeable difference was the total. While the S1 scenario uses only twice the estimated best-fit volume of 5.6 Mm³, the volume for S3 is six times larger, and the one for S3 with the extended sectors is 10 times larger. In contrast, the volume estimated from field observations in Cole et al. (2023) at 17 Mm³ is pretty much the centre of our volume range for the rapid simulations (10 < V < 30 Mm³). Particular care must be taken regarding the total volume estimation of this best-fit simulation, owing to the simplified source conditions used (single-pulse event) that could explain the discrepancies observed between the total simulated flow volume v. the total deposit volume as estimated from field observations. This reinforces the idea that the total volume of a PDC is the most challenging input parameter to constrain. Very few data can help estimate a volume, but they are not always reliable. For instance, the volume of a lava dome can be used as a volume proxy for block-and-ash flows (Kelfoun et al. 2017; Ogburn and Calder 2017; Gueugneau et al. 2019), but there is no evidence for a total or partial dome collapse, or if more material will be incorporated in the PDCs by a subsequent explosion. This is why we decided to use an average volume for common small-volume PDCs from FlowDat (Ogburn 2012), 10 Mm³, and to double and triple it to investigate its impacts on PDC extent during the syn-crisis phase. Such volume overestimation allows us to encompass all potential areas that might be at risk of PDC invasion. Thanks to that, no false negative, i.e. area inundated by the real PDCs left uncovered in our simulations, was found between the syn-crisis cumulative maps and the field deposit map (Figs 1c, 3–5) and the volume estimate from Cole et al. (2023) demonstrates that our syn-crisis volumes were very realistic.

The comparison between the syn-crisis and post-crisis simulations also highlights the successful choices made when selecting inputs for all simulation scenarios, and the ability of the code to accurately model small-volume PDCs. This is also the case for most of VolcFlow's input parameters, such as the particle properties (i.e. particle size and density, gas properties, etc.) that remained mostly unchanged between the two sets of simulations. It is also true for the source properties (position and shape), which were correctly chosen since the first simulation owing to the known position of the new crater, even though the shape had to be slightly refined for the post-crisis simulation. Positioning the source along the crater rim at the beginning of each inundated valley was found to be the best configuration and fits well with a fountain or a column collapsing inside a specific sector of the crater, for a short time period.

To summarize, this exercise demonstrates that VolcFlow has the potential to accurately model such PDCs during a crisis when the correct input and source parameters are selected. Data from this eruption are a new addition to VolcFlow input parameters database that can be used for future crises to reduce the preparation time or help to define ranges of value to explore in the case of probabilistic hazard assessment. Most of the work should now focus on better constraining input parameters as soon as an eruption begins, which relies on (1) discretizing each PDC pulse in time and space to facilitate the determination of the most at-risk zones; and (2) estimating their parameters simultaneously (volume, inundated area, mass flux) to avoid the use of simplified source conditions. Continuous visual monitoring of volcanic edifices, such as at Merapi for example (Kelfoun et al. 2021), or continuous satellite image acquisition as discussed in the next section, seem to be promising methods.

Lack of knowledge of input parameters is the main limitation of a deterministic approach, especially in the case of rapid modelling. The main source of data is usually visual observations and satellite imagery. However, in the case of the April 2021 La Soufrière eruption, visual observations and satellite imagery were extremely limited owing to the close spacing of explosions and intense ash fallout from the early parts of explosive activity. Radar imagery appears to be an effective solution to tackle visibility issues since radar wavelengths are able to pass through either ash or clouds. Here, radar images provided important data: (1) the position of the crater from radar images (Fig. 2a); (2) the evolution of the topography (Fig. 2b), which could lead to new DEMs; and (3) the tracking of PDC deposits with radar backscatter imagery (Fig. 2c), to ensure a correspondence between the ongoing simulations and PDCs forming during the eruption.

Using the backscatter signal from satellite radar images acquired during an eruption, like the ones used here from TerraSAR-X (Fig. 2c), has the potential to provide useful information about the eruptive processes happening both at the source and on the flanks of the volcano. However, backscatter changes are complex and challenging to interpret: explosive deposits produce different signals depending on pre-existing ground cover, radar parameters and eruption characteristics, their thickness varies over several orders of magnitude and they tend to be rapidly remobilized and eroded (e.g. Wadge et al. 2011; Solikhin et al. 2015; Dualeh et al. 2021). The characterization of subtle backscatter signals associated with explosive eruptions are best observed when using a combination of (1) high temporal- and spatial-resolution backscatter imagery, (2) radiometric terrain calibration, (3) speckle correction and (4) consideration of pre-existing scattering properties. However, the interpretation of synthetic aperture radar backscatter for volcanology is challenging because there is no simple relationship between the magnitude or sign of backscatter change and the physical properties of fresh volcanic deposits. For volcanic deposits, the surface roughness generally contributes the most change to the radar backscatter signal recorded (Wadge et al. 2011; Solikhin et al. 2015). Changes in backscatter also depend strongly on the scattering properties of the previous surface cover, resulting in complex change patterns on the vegetated volcano flanks (e.g. southwestern sector; Fig. 2c). For example, where a PDC removes vegetation, the ground becomes smoother on the scale of the X-band radar wavelength (3.1 cm) and the contribution of volumetric scattering is removed, resulting in a decrease in backscatter (blue areas in drainages and interfluves of the southwestern flank; Fig. 2c). Over densely vegetated areas of the volcano flanks, tephra fallout can also cause a decrease in backscatter, but appears to be more limited possibly owing to variations in ash thickness, ash and ground moisture content and pre-eruption land cover. Finally, using this approach does not allow us to capture (1) low-magnitude backscatter changes that were only slightly larger than background variations (i.e. lower flanks not affected by PDCs), (2) areas where backscatter changes from different deposits overlap (i.e. summit area) or (3) areas where the valley-confined PDC deposits are emplaced inside relatively small and narrow channels (i.e. upper SW valleys). Using a dual-band approach (i.e. NASA ASAR-ISRO L + S band data), or combining multiple polarimetric channels, could allow for some of these features to be investigated and mapped with backscatter, since explosive deposits of different grain sizes and sorting will produce different backscatter signals at specific band wavelengths (Rogic et al. 2023). This will be especially the case for different PDC pulses deposited in the same zone but sufficiently well separated in time (owing to the repeat time of the satellite used, usually a few hours at least), for which isolating their own backscatter signal could help discretize each PDC pulse in near-real time.

Most eruptions at St Vincent have generated PDCs, and it is worth comparing the extent of those with the simulations undertaken in this study. The 1979 eruption generated PDCs of quite similar extent to the April 2021 eruption (VEI 3 v. VEI 4; see Cole et al. 2019, 2023), reaching the sea at the Larikai and Roseau valleys and generating a limited ash-cloud surge component in proximal regions (Shepherd et al. 1979). The 7 May 1902 eruption, however, generated much more extensive PDCs (Anderson and Flett 1903). Thick concentrated PDCs were contained within the numerous valleys including the Wallibou, Dry Wallibou, Roseau and Rabacca; however, the dilute ash-cloud surge component was extensive (Fig. 8). Contemporary 1902 reports (Anderson and Flett 1903) describe people in boats a few miles offshore to the west of the volcano being engulfed by the ash-cloud surge. On the southeastern side, thick concentrated PDCs filled the Rabacca valley, changing the geomorphology, although PDCs did not reach the sea in this region. The ash-cloud surge inundated the lower part of Richmond Vale, which is separated from the Wallibou valley draining the volcano by a ridge (Bunkers hill). Thus, it appears the ash-cloud surge swept around to the south on exiting the mouth of the Wallibou valley and travelling over the sea. The inundation area of the ash-cloud surges formed in 1902 are quite similar in some ways to the ash-cloud surges simulated for the ‘worst case scenarios’ in Figure 5. The area inundated to the SW in 1902 probably fits best with the S2 simulation, even though only the concentrated part of 1902 PDCs reached the sea. The fit is less similar in the worst-case scenarios to the SE as PDCs in 1902 did not reach the sea at the mouth of the Rabacca. Nevertheless, ash-cloud surge inundation is most similar to the S1 simulation. The volume estimate of Hay (1959) of about 9.9 Mm³ (3.5 × 10⁸ cubic feet DRE, converted in bulk volume) is almost identical to our S1 simulation, and supports again that the worst-case scenarios defined here (Fig. 5) are certainly possible. Nevertheless, the same care regarding the single pulse simplification must be considered with the 1902 eruption as for the 2021 eruption.

A crucial but challenging step in such a rapid scenario-based project was to determine in the smallest amount of time possible the best balance between the DEM resolution and simulation time. The highest resolution DEM possible is always the goal for such simulations, to ensure precise and reliable results; however, this can have a dramatic effect on the calculation time. For instance, with the same computational domain as that in Figures 4 or 6, it takes c. 3 h to run a two-phase VolcFlow simulation on a regular desktop computer (CPU i7-4400 K) over the 8 m resolution DEM, but this rises to 1–2 days using the same DEM at 2 m resolution. However, it took us the entire day of 10 April to come to the conclusion that 8 m was the most suitable DEM resolution. This duration is not negligible, but pre-eruption planning could reduce such time loss. Determining the highest risk sector and estimating the most suitable DEM resolution in respect to the size of the computational domain before an eruption could drastically reduce the preparation time and make the code immediately operational for rapid hazard assessment. It would also give immediate access to the DEM, made ready-to-use, and reduce the time spent by finding and getting access to the most recent DEM with the highest resolution possible.

Another challenge concerns the uncertainty quantification that cannot be determined during the syn-crisis phase. Non-negligible epistemic and aleatory uncertainties are associated with the poorly constrained input parameters and the natural variability of the eruptive system. Quantification of these uncertainties was only possible in post-crisis simulations with Bayesian metrics using well-constrained field data, implying a comparison between a cumulative PDC deposit field map resulting from the emplacement of multiple PDC pulses over 2 weeks with a map obtained from the simulation of a single PDC pulse, which obviously contains some uncertainties. Without a solution for such uncertainty quantification, further model validation of other well-studied PDC can help decrease uncertainties by reducing the range of the code's input parameters.

Overall, the results of this work demonstrate that such a rapid scenario-based hazard assessment study is an excellent complement to more conventional, long-term probabilistic strategy and brings important insights that could not be obtained otherwise. First, this study used, for the first time, a two-phase depth-averaged model during an eruptive crisis. This type of model differs from empirical/kinetic models typically used for hazard assessment, such as ECM (the Energy Cone Model; Malin and Sheridan 1982), LAHARZ (Schilling 1998), the Box model (Esposti Ongaro et al. 2016) or more recently ECMapPro (Aravena et al. 2020). Since the computation of flow kinematics with our model is physically based, the complex interaction between PDCs and the topography, as recent studies have highlighted (Charbonnier and Gertisser 2008; Lube et al. 2011; Ogburn et al. 2014; Kubo Hutchison and Dufek 2021), can be accurately reproduced, in contrast to kinetics model that tend to underestimate the effects of topography (Ogburn and Calder 2017; Sandri et al. 2018; Aravena et al. 2020). In addition, the two-phase version allows for the modelling of the two regimes of PDCs, i.e. the concentrated flow and the ash-cloud surge. Each regime can be represented either separately as cumulative maps (Figs 4–6) or together as a post-simulation map (Fig. 7), which brings a new dimension for PDC hazard assessment as demonstrated at El Misti (Charbonnier et al. 2020). As a consequence, we are able to infer which type of PDC regime a specific area will be mostly affected by. We can also test realistic scenarios based on eruptive processes, i.e. fountain or dome collapse, leading to a more specific PDC forecast, or estimate the volume of a PDC during an ongoing eruption, for example. This approach is suitable for future volcanic crises, allowing resources to be focused on a syn-crisis hazard assessment, as suggested by Newhall and Pallister (2015) or Rouwet et al. (2019), aiming at testing specific scenarios, or ensuring real-time PDC forecasting, as well as complementing probabilistic hazard mapping.

We have presented a two-step modelling study of PDCs generated during the April 2021 La Soufrière, St Vincent eruption with the two-phase version of VolcFlow, initiated to support the response team. Preliminary simulations run the day of the first explosion on 10 April highlighted that the southwestern and the southeastern flanks of La Soufrière were the areas with the highest potential for PDC invasion. The first modelling step conducted during the eruptive phase, focusing on the evaluation of PDC invasion in these two sectors, was based on the simulation of three PDC scenarios: 10, 20 and 30 Mm³ PDCs. Cumulative maps summarizing simulation results confirmed the preliminary simulations and showed that: (1) the southwestern sector is at most risk of PDC invasion; (2) PDCs are well channelized by La Soufrière deep valleys and can easily reach the sea on the west coast, but also on the east coast for the most voluminous scenarios; (3) Chateaubelair and Georgetown can be reached by a PDC but only in the case of an extremely voluminous PDC (>30 Mm³), not likely to occur with the current activity of the volcano at the time; and (4) while the PDC hazards are satisfactorily estimated in the ‘very high risk’ zone of the volcanic risk map of St Vincent, the Chateaubelair sector must be closely monitored. In the second modelling step conducted months after the crisis, the newly acquired data from the eruption, such as a PDC deposit map, were used for a new set of simulations aiming at best reproducing the observed PDCs. Input parameters used for the best-fit simulation are then compared with those of the syn-crisis simulation. The results showed that ash-cloud surges in syn-crisis simulations were overestimated (two order of magnitude for the surge production coefficient), and that the volume was also too high (5.6 Mm³ for the best-fit simulation, half of the S1 scenario). They also highlighted that our source positioned along the crater rim was accurate, and that most of the input parameters of VolcFlow were already correctly estimated. Finally, the single PDC pulse modelling approach employed here did not fully reproduce the characteristics of the real eruptive sequence comprising multiple PDC pulses, but it has proven to be a robust substitute in the case of a lack of information about the source conditions as it satisfactorily approximated the cumulated deposit map at first order.

This study displays the first example of rapid scenario-based assessment with a two-phase depth-averaged model, and discusses for the first time the strengths and weakness of the choice of poorly to unknown input parameters. Important lessons were learned on: (1) how to best estimate and represent the PDC invasion potential (cumulative maps); (2) how crucial the DEM resolution selection is; (3) how radar images can provide essential syn-crisis data; and (4) how VolcFlow's simulation results can be improved when input parameters are well constrained. The new set of input parameters used for the best-fit simulation is a new addition to a database constituted of other PDC modelling projects with this code. The different sets of input parameters help to define a range of values that is increasingly precise and can be used for future similar PDC hazard assessment projects at other volcanoes, prior to or during an eruption.

The authors would like to thank all personnel involved in the response team that provided feedback on our cumulative maps, as well as UWI-SRC and the Montserrat Volcano Observatory for providing us with the DEM of St Vincent Island. We also thank R. Richardson for his editorial work.

The authors declare no competing interests.

VG: conceptualization (lead), data curation (lead), formal analysis (lead), investigation (lead), methodology (lead), project administration (equal), software (lead), writing – original draft (lead); SC: conceptualization (supporting), formal analysis (supporting), funding acquisition (lead), supervision (supporting), validation (supporting), writing – review & editing (supporting); VLM: conceptualization (supporting), data curation (supporting), formal analysis (supporting), investigation (supporting), project administration (supporting), validation (supporting), visualization (supporting), writing – original draft (supporting), writing – review & editing (supporting); PC: investigation (supporting), resources (supporting), visualization (supporting), writing – original draft (supporting), writing – review & editing (supporting); RG: resources (supporting), writing – review & editing (supporting); EWD: resources (supporting), writing – review & editing (supporting).

This work was supported by the National Science Foundation CAREER grant no. 1751905. Satellite acquisitions were supported by public funds received in the framework of GEOSUD, a project (ANR-10-EQPX-20) of the program ‘Investissements d'Avenir’ managed by the French National Research Agency.

VolcFlow Two-phase is freely available at https://lmv.uca.fr/volcflow/. Simulation outputs, input files and movies are available upon request from the corresponding author. Satellite images were provided by R. Grandin (Capella Space) and the TerraSAR-X imagery (acquired by DLR) and can be requested through the CEOS Volcano Demonstrator Program, by contacting either the demonstrator lead or the WGDisasters secretariat (https://ceos.org/ourwork/workinggroups/disasters/). Access to Pleiades data was granted through the DINAMIS program (https://dinamis.teledetection.fr/ via project ID 2021-055-Sci (PI: Raphael Grandin, IPGP, Grandin and Delorme 2021), and via the CIEST2 program for ForM@Ter (https://www.poleterresolide.fr/ciest-2-nouvelle-generation-2/). Calculation of the Pleiades DEM used the S-CAPAD cluster of IPGP.