Modelling pyroclastic density currents of the April 2021 La Soufriere St. Vincent eruption: from rapid invasion maps to field-constrained numerical simulations

.


Introduction
After 41 years of quiescence, two 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 April 9 th , 2021, and entered in a new eruptive phase that lasted two weeks.The response team, assembled at the time of the first extrusive activity, successfully forecasted the transition to explosive activity, enabling evacuation of areas at the highest risk around La Soufrière on April 8 th .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 phenomenon that could result from the eruption of La Soufrière and the destruction of the newly formed lava dome.Among them, pyroclastic density currents (PDC) 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 have already occurred during the 1902-1903(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-1903(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 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 modeling 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): (i) 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, (ii) 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; A C C E P T E D M A N U S C R I P T Downloaded from https://www.lyellcollection.orgby University of Plymouth, Charles Seale-Hayne Library on Jun 05, 2023 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 guaranties a more complete hazard forecasting at many levels (Newhall and Pallister, 2015;Rouwet et al., 2019).Regarding La Soufrière of St. Vincent, probabilistic hazard assessment has only been recently investigated by 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 (i) assessing PDC hazards during and after the April 2021 eruption, and (ii) better constraining input parameters of numerical simulations to improve rapid syn-crisis hazard assessments generally.The study is therefore organized in two steps: (i) 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 (ii) a post-crisis phase during which syn-crisis maps were reanalyzed 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 two weeks after the onset of explosive activity on 9 th April.We chose this code for: (i) its capacity to simulate the extent of both dilute and concentrated regimes of PDCs during a single simulation, (ii) its rapidity, i.e., simulation run < 5h, (iii) the simplicity to generate invasion maps of concentrated and/or dilute PDCs from simulation outputs.Cumulative PDC invasion maps are 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, difference in data processing, capability of the two-phase version of VolcFlow for rapid hazard assessment, and the pros and cons of scenario-based approach for PDC hazard assessment, are discussed in terms of lessons learned for future eruptions.

La Soufrière volcano
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 ~1500 fatalities and a considerable economic impact.The last eruption in 1979 produced a series of thirteen explosions and resumed with a lava dome growth at the center (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 (FIRMS) detected a thermal anomaly inside the summit crater, which was confirmed as new dome growth by visual observations on 29 December 2020 (Miller et al., 2022).This effusive activity continued for a period of three 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 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 (Appendix, this volume).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 preexisting 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 (Miller et al., 2022).Cole et al., 2023;this issue).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.

Pyroclastic density currents of the April 2021 explosions
A precise timeline of all PDC events could not be created due to the lack of visual observations.For the first 48 hours of the eruption, copious ash in the atmosphere as a result of closely spaced explosions caused darkness, making visual observations difficult or impossible.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 of fine lapilli, typically poor in fine ash.These deposits are relatively continuous but do show some minor lateral thickness variations.

Method: the numerical model and the DEMs used
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 modeled using depthaveraged 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 due to 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 depthaveraged approach, the plastic rheological law is used, involving a constant retarding stress T , which have been shown to successfully reproduce various features of such currents and their deposits (Kelfoun et al., 2009(Kelfoun et al., , 2017;;Kelfoun, 2011;Charbonnier and Gertisser, 2012;Ogburn and Calder, 2017;Gueugneau et al., 2019Gueugneau et al., , 2020)).The ash-cloud surge is simulated as a turbulent continuum that loses momentum due 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 th April -R.Grandin, IPGP), downgraded to 8 m resolution to allow faster simulation runs.Due 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 The University of the West Indies Seismic Research Centre (UWI-SRC) in 2005 (Robertson, 2005).The merged DEM was then resized to 8 m cells for consistency with the first simulations.

Preliminary simulation, 10 th April 2021
We ran preliminary simulations on April 10 th 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 th April (Fig 2A), the position of the new crater was located, and test simulations could be run with a source centered 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-or 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., 2019Gueugneau et al., , 2020)), the value 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 3 ) 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 due to major topographic barriers preventing their northwards propagation.
These conclusions and our maps in Fig. 3 were shared with UWI-SRC, and two different sectors were selected for further investigation: (i) the southwestern flank, with the highest probability of inundation based on the preliminary simulation, (ii) 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 Fig. 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.

Rapid scenario-based approach: procedure and parameters
To generate an inundation map, a scenario was defined, and a series of simulations was run for which the value 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 parameters were estimated and calibrated from values used in previous studies (Table 1).The last two parameters of the table, a.k.a.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 twelve other parameters fixed in all simulations.Three different volumes were explored to define three PDC-type scenarios: S1, with a volume of 10 Mm 3 , typical of small-volume PDCs (FlowDat database, Ogburn, 2012), similar to the June 25 th , 1997 PDCs at Soufriere Hills, Montserrat (Loughlin et al., 2002); (ii) S2, with a volume of 20 Mm 3 , twice the targeted volume ; (iii) S3, with a volume of 30 Mm 3 , three times the targeted volume, similar to the May 8 th , 1902 PDCs at Mount Pelée, Martinique (Lacroix, 1904;Gueugneau et al., 2020).A supply duration of 75 s was chosen to ensure a volumetric rate at the source in the range of 0.1 to 0.5 Mm 3 s -1 , typical of small-volume PDCs (Loughlin et al., 2002;Gueugneau et al., 2020;Roche et al., 2021).The position of the source was determined by our knowledge of the site of the April 9 th explosion, and the position of the lava dome since December 2020.Note that PDCs are modeled here as a single-pulse collapse event in each scenario, rather than multiple pulses resulting from multiple collapsing events at the source, as it seems to have been the case during this eruption.Simplifying such complex process into one single pulse can affect the reliability of our simulations and maps, but it also reduced the number of source parameters (and their uncertainties) used to simulate such PDC events.

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

Southeastern sector, (12 th April simulation)
To evaluate the PDC invasion potential on the southeastern sector of La Soufrière, we set an arcuateshape source along the southeastern quarter of the crater rim (see Fig. -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 3 ) that produces an extensive ash-cloud surge.Otherwise, ashcloud surges seem restricted inshore.

Southwestern sector, April 14 th
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.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 toward the sea.Therefore, Chateaubelair seems out of reach, and only populations around the Wallibou valley's mouth are threatened by the passage of a PDC, and only for very voluminous scenarios (S2 and S3).
Here again, the effect of an increasing volumetric rate is mostly restricted to an increase of 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 to conclude 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.

Extended sectors, Chateaubelair and Georgetown
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 km 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's sector (left map in Fig. 6), the source was  (Edmonds and Herd, 2005), at Krakatoa in Indonesia (Carey et al., 1996), or also numerically at Mt 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 Fig. 4 to restrict the focus on the Rabacca drainage system only.Similarly, only voluminous scenarios (i.e., >30 Mm 3 ) are able to produce a surge that inundates Georgetown, with an enhanced lateral spreading toward the sea.We note however that the Mt Brisbane hills are smoother than the Richmond hills from the southwest sector, and oriented in a NE-SW direction rather that E-W, 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 3 ), which is not compatible with a Vulcanian or violent Vulcanian type PDCs, but more a blast-like PDC (Mt Pelée in 1902, Gueugneau et al., 2020;Montserrat in 1997, Sparks et al., 2002) 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 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 Figs 4 and 5, can be potential emerging points for PDCs.Our maps allowed us to: (i) identify sectors at most risk, (ii) identify the worst case scenarios leading to an invasion of Chateaubelair and Georgetown, (iii) 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, weather 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.

Procedure and parameters
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., 2022, this volume).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 (Asim) and the area of the PDC deposit in the field (Aobs): 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 (RJ) uses TP and the union of areas inundated by observed and simulated flows: • the model sensitivity (RMS) uses TP and the area inundated by simulated flows: • the model precision (RMP) uses TP and the area inundated by observed flows: A trial-and-error method was adopted, and input parameters were varied, one at a time, until we obtained the highest metric values possible that best reproduced the extent of the real PDCs.The value of each VolcFlow input parameter is reported in Table 2, in comparison with the values we selected for the syn-crisis simulations.Results of the validation metrics calculation for the best-fit simulation, expressed as percentage of match, are also reported at the bottom of the Table 2.Note that the field deposit map is made of an accumulation of several PDCs pulses emplaced during the two weeks of the crisis, and while these new field data provided important insights about the overall PDC sequence (Cole et al. 2023), many uncertainties remain about the deposit characteristics of each individual pulse .Therefore, similar source conditions as those used for the syn-crisis simulations are applied here.

Best fit simulation
The result of the best-fit simulation is presented in Fig. 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 due to the low resolution of the DEM (8 m), resulting in relatively low values for the 3 metrics in the Table 2.In addition, the longer flow runout in upper Wallibou valley did not matched the actual PDC.The simulated ash-cloud surge, its extent is quite close to the observed one in Fig. 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 due 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-shape 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 disk 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 the 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 3 , half of the volume of S1, the smallest-volume scenario.For finer adjustments, the constant retarding stress (CRS) of the concentrated flow was increased to obtain shorter runouts, in accordance with previous studies showing that the CRS is inversely proportional to the flow volume (Kelfoun, 2011;Charbonnier et al., 2020).

Inter-comparison between rapid maps and the new simulations
The 'post crisis' best fit simulation map presented in Fig. 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 that the total.While the S1 scenario uses only two times the estimated best-fit volume of 5.6 Mm 3 , the volume for S3 is six times larger, and the one for S3 with the extended sectors is ten times larger.In contrast, the volume estimated from field observations in Cole et al. (2023) at 17 Mm 3 is pretty much the center of our volume range for the rapid simulations (10<V<30 Mm 3 ).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 versus 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 following explosion.This is why we decided to use an average volume for common small-volume PDCs from FlowDat (Ogburn, 2012), 10 Mm 3 , 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, were found between the syn-crisis cumulative maps and the field deposit map (Fig. 1C, 3, 4 and 5).And the volume estimate from Cole et al. (2023) demonstrate 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…) that remained mostly unchanged between the two sets of simulations.It is also true for the source properties (position and shape), that 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 crisis to reduce the preparation time or help to define ranges of value to explore in 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 (i) discretizing each PDC pulse in time and space to facilitate the determination of the most at risk zones; and (ii) 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 paragraph, seem to be promising methods.

Satellite images: the importance of radar imagery
Lack of knowledge of input parameters is the main limitation of a deterministic approach, especially in the case of rapid modeling.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 due 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: (i) the position of the crater from radar images (Fig. 2A); (ii) the evolution of the topography (Fig. 2B), which could lead to new DEMs; (iii) 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 SAR 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 appear to be more limited possibly due 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 (a) low magnitude backscatter changes that were only slightly larger than background variations (i.e., lower flanks not affected by PDCs), (b) areas where backscatter changes from different deposits overlap (i.e., summit area) or (c) areas where the valley-confined PDC deposits are emplaced inside relatively small and narrow channels (i.e., upper southwest 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., under review).It 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 discretizing each PDC pulse in near-real time.(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.

Comparison of Simulations and
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 Fig. 6.The area inundated to the southwest in 1902 probably fits best with the S2 simulation.The fit is less similar in the worst-case scenarios to the southeast 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 ( 1956) of about 9.9 Mm 3 (3.5 10 8 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. 6) are certainly possible.Nevertheless, the same care regarding the single pulse simplification must be considered with the 1902 eruption as for the 2021 eruption.

Lessons learned for rapid scenario-based PDC hazard assessment
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 the one in Fig. 4 or 5, it takes ~3 hr to run a two-phase VolcFlow simulation on a regular desktop computer (CPU i7-4400K) 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 th 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 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 twophase depth-averaged model during an eruptive crisis.This type of model differs from empirical/kinetic models typically used for hazard assessment, such as ECM (Energy Cone Model; Malin and Sheridan 1982), LAHARZ (Schilling, 1998), 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 Hutchinson 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 modeling 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 (Fig. 4, 5 and 6), or together as a post-simulation map (Fig. 7), which bring 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 to 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 a real-time PDC forecasting, as well as complementing probabilistic hazard mapping.

Conclusion
We evaluation of PDC invasion in these two sectors, was based on the simulation of three PDC scenarios: 10 Mm 3 , 20 Mm 3 , and 30 Mm 3 PDC.Cumulative maps summarizing simulation results confirmed the preliminary simulations and showed that (i) the southwestern sector is at most risk of PDC invasion, (ii) PDC 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, (iii) Chateaubelair and Georgetown can be reach by a PDC but only in case of extremely voluminous PDC (> 30 Mm 3 ), not likely to occur with the current activity of the volcano at the time, and (iv) 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 modeling 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 bestfit simulation are then compared with those of the syn-crisis simulation.Results showed that ashcloud 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 3 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 modeling approach employed here did not fully reproduce the characteristics of the real eruptive sequence comprised of multiple PDC pulses, but it has proven to be a robust substitute in case of 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  showing the new crater formed after the April 9th first explosion, and its morphological evolution after two weeks of activity.Note the thick tephra deposit that piles up along the northern side of the new vent.c) TerraSAR X radar backscattering image difference between April 1st and April 12th, geocoded and multilooked.Red color represents an increase in the backscattering, yellow a neutral value, and blue a decrease.The white line materializes the ash-cloud surge extent as observed in the field (see Figure 1).The blue disk represents the area where ash-cloud surge is supplied, and the yellow croissants crescents represent the area where both concentrated flow and surge are supplied.For ease of comparison, the deposit map is added in the top left corner.
://www.lyellcollection.orgby University of Plymouth, Charles Seale-Hayne Library on Jun 05, 2023 ://www.lyellcollection.orgby University of Plymouth, Charles Seale-Hayne Library on Jun 05, 2023 PDCs were associated with the explosive activity and emplaced mainly onto the southwestern flanks of the volcano.The PDC deposit map created after field campaigns in early 2022 is presented in Fig 1C (see also Satellite imagery (radar backscattering, Fig 2C) indicated that the first PDCs had been emplaced by 11 th 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 th 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 (~1 km) down valleys to the southeast (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 3 to 5 PDC units of 1 to 4 m thick, giving a total deposit thickness of 20 to 30 m. Dilute PDCs, or ash-cloud surge, were extensive in proximal regions (Fig 1C, light red area), ://www.lyellcollection.orgby University of Plymouth, Charles Seale-Hayne Library on Jun 05, 2023 ://www.lyellcollection.orgby University of Plymouth, Charles Seale-Hayne Library on Jun 05, 2023 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 toward 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 midpoint, 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 toward the sea as Mount Brisbane hills block their southward propagation.Differences in flow extent between the three scenarios are important with an increase of the inundated area by a factor of two between S1 and S3.Although the modeled 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: from https://www.lyellcollection.orgby University of Plymouth, Charles Seale-Hayne Library on Jun 05, 2023 ://www.lyellcollection.orgby University of Plymouth, Charles Seale-Hayne Library on Jun 05, 2023 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 3 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 toward the coast as the relief decreases and reaches the northern end of Chateaubelair.We then increased the volume further to 50 Mm 3 for S3, which causes 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, modeled 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 ://www.lyellcollection.orgby University of Plymouth, Charles Seale-Hayne Library on Jun 05, 2023 ://www.lyellcollection.orgby University of Plymouth, Charles Seale-Hayne Library on Jun 05, 2023 ://www.lyellcollection.orgby University of Plymouth, Charles Seale-Hayne Library on Jun 05, 2023 ://www.lyellcollection.orgby University of Plymouth, Charles Seale-Hayne Library on Jun 05, 2023 PDCs formed in past eruptions at St. Vincent 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 4 vs VEI 5; 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 ://www.lyellcollection.orgby University of Plymouth, Charles Seale-Hayne Library on Jun 05, 2023 ://www.lyellcollection.orgby University of Plymouth, Charles Seale-Hayne Library on Jun 05, 2023 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 wellconstrained field data, implying a comparison between a cumulative PDC deposit field map resulting from the emplacement of multiple PDC pulses over two 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.
have presented a two-step modeling study of PDCs generated during the April 2021 La Soufrière of 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 April 10 th highlighted that the southwestern and the southeastern flanks of La Soufrière were the areas with the highest potential of PDC invasion.The first modeling step conducted during the eruptive phase, focusing on the ://www.lyellcollection.orgby University of Plymouth, Charles Seale-Hayne Library on Jun 05, 2023

Figure Captions Figure 1 :
Figure Captions

Figure 2 :
Figure 2: a & b) Radar amplitude images from Capella Space taken on April 10th and April 23thshowing the new crater formed after the April 9th first explosion, and its morphological evolution after two weeks of activity.Note the thick tephra deposit that piles up along the northern side of the new vent.c) TerraSAR X radar backscattering image difference between April 1st and April 12th, geocoded and multilooked.Red color represents an increase in the backscattering, yellow a neutral value, and blue a decrease.The white line materializes the ash-cloud surge extent as observed in the field (see Figure1).

Figure 3 :
Figure3: Results of the preliminary simulation ran on April 11th: A) final state of the simulation, where the position of the source is represented by the orange disk, the concentrated deposits in red/purple, and the ash-cloud surge deposits in gray/green.B) to E) snap shots of the simulation in 3D at22, 68, 128, and 192  s showing the emplacement dynamics of each layer (concentrated flow in red, ash-cloud surge in blue until it lifts off and only displays its deposits in green).

Figure 4 :
Figure 4: Cumulative PDC invasion maps (bottom) for the southeastern corner of La Soufrière volcano, created and distributed to UWI-SRC on April 12th.Cumulative maps are generated by superimposing syn-crisis simulation results (simulations S1, S2 and S3, top row) of either the concentrated flow (bottom left) or the ash-cloud surge (bottom right).The source used to generate the simulated flows is materialized by the black arc along the crater rim.

Figure 5 :
Figure 5: Cumulative PDC invasion maps (bottom) for the southwestern corner of the Soufrière volcano, created and distributed to UWI-SRC on April 13th.Cumulative maps are generated by superimposing syn-crisis simulation results (simulations S1, S2 and S3, top raw) of either the concentrated flow (bottom left) or the ash-cloud surge (bottom right).The source used to generate the simulated flows is materialized by the black arc along the crater rim.

Figure 6 :
Figure 6: Cumulative PDC invasion maps for Chateaubelair (A) and Georgetown sectors (B) using the extended DEM, created and distributed to UWI-SRC on April 13th.Maps are based on the ashcloud surge extent only.To investigate the worst case scenarios, the initial volume for the 3 scenarios were set to 10, 30 and 50 Mm3.

Figure 7 :
Figure 7: Results of VolcFlow's best fit post-crisis simulation, compared to the field deposit map shown in Fig. 1.Concentrated deposits are represented by the purple to red colormap, and surge deposits by the gray to green colormap, with the associated thicknesses shown under the colorbars in the legend box.The black line materializes indicates the extent of the surge as mapped in the field.The blue disk represents the area where ash-cloud surge is supplied, and the yellow croissants crescents represent the area where both concentrated flow and surge are supplied.For ease of comparison, the deposit map is added in the top left corner.