Volcanic crisis management supported by near real-time lava ﬂow hazard assessment at Piton de la Fournaise, La Réunion

Since 1979, Piton de la Fournaise (La Réunion) has erupted on average two times per year, with 95 % of these eruptions occurring within an uninhabited caldera. However, lava ﬂows have occasionally impacted populated regions on the island, as in 1977 and 1986. Since 2014, an integrated satellite data–driven multinational response to effusive crises has been developed to rapidly assess lava inundation area and ﬂow runout distance. In 2018, this protocol was implemented as a standalone software to provide a lava ﬂow hazard map showing the probability of ﬂow coverage and runouts as a function of discharge rate. Since 2019, the produced short-term hazard map is shared with local civil protection in the ﬁrst few hours following the start of an eruption to aid in mitigation actions. Multiple exchanges between scientists, the observatory, and civil protection has improved the delivered hazard map, ensuring a common understanding, a product which is of use and usable, and helping to build effective mitigation strategies at Piton de la Fournaise. In this work we illustrate this effective near real-time protocol with case studies and document how the produced short-term hazard map has been tailored to meet the needs of civil protection.


INTRODUCTION
Lava flows are one of the most frequent volcanic hazards in the world, yet they are considered to rarely pose a risk to people [Witham 2005;Harris 2015]. Indeed, at most sites, lava flows are slow enough to allow populations to be evacuated [e.g. Weisel and Stapleton 1992;Duncan et al. 1996;Chester et al. 1999]. However, examples where lava flows have caused fatalities do exist [Blong 1984], such as during the flank eruptions of Nyiragongo in the Democratic Republic of Congo in 1977Congo in , 2002Congo in , and 2021. Lavas associated with these Nyiragongo eruptions were extremely fluid [Giordano et al. 2007;Morrison et al. 2020] and thus advanced extremely rapidly (at tens of km h −1 [Tazieff 1977]). As a result, lava flows reached the city of Goma within a few hours, causing dozens of fatalities either directly (burns, asphyxia) or indirectly (e.g. explosion of a gas station or car accidents during evacuation) [Harris 2015].
However, by far the greatest problem facing a community in the path of a lava flow is the destruction of structures by burning, flooding, and/or burial [Blong 1984]. Any unmovable structure, and its contents, will be destroyed including buildings, businesses (factories, garages, shops), essential community facilities (schools, churches, sports centers, hospitals), infrastructure, utilities, and agriculture [e.g. Finch and Macdonald 1950;Luhr and Simkin (eds) 1993;Duncan et al. 1996]. Destruction of all of these essential elements of a community have lasting repercussions on social fabric, mental health, and local economies [Harris 2015]. For example, after the 2002 Nyiragongo eruption, about 120,000 people were left homeless, and 15 % of the city and a third of the airport runway were completely buried by lava [Tedesco et al. 2007]. Following the 2021 Nyiragongo eruption, 3600 houses were de-stroyed and more than 400,000 people were affected either by water supply difficulties or displacement (UNOCHA report of 5 June, 2021 * ). In June 2018, the effusive eruption of Kīlauea in the Puna region of Hawai ' i destroyed more than 700 structures and caused economic losses estimated at more than 800 million US dollars † [Meredith et al. 2022]. Finally, the eruption of La Palma (Canaries), which began in September 2021 and ended in January 2022, buried ca. 12 km 2 of land including agricultural fields and villages, destroying more than 1600 houses ‡ . As a result, more than 7000 people were displaced § .
To better assess lava flow hazard, run real-time appraisals of potential lava flow inundation areas and, thus, respond to an effusive crisis in an effective manner, scientists and civil protection agents must work together closely to test and validate operational tools before an event occurs. With an average of two, mainly effusive, eruptions per year over the last 40 years, Piton de la Fournaise volcano (La Réunion, France) is a perfect laboratory to develop and test effusive crisis response protocols . Piton de la Fournaise is an active basaltic hotspot volcano on La Réunion Island, a French overseas department located in the Western Indian Ocean (Figure 1A). Since 2014, scientists have been working on an integrated satellite data-driven response for effusive crises at Piton de la Fournaise, with the aim of assessing, in near realtime, potential lava inundation areas and flow runout distances . This protocol has now been implemented as a standalone software package that can be tied to any Geographic Information System (GIS) allowing production of a short-term hazard maps within the first minutes to hours after the start of an eruption ]. Maps are a common way of communicating comprehensive hazard information allowing the spatial extent of the hazard, as well as information as to key infrastructure and land type at risk, to be quickly assessed by emergency managers [Pareschi et al. 2000; Thompson et al. 2015].
In the present study, we detail this near real-time protocol, package, and product. We first provide background on mitigating the risks associated with lava flows at effusive centers in general, and then focus on crisis management at Piton de la Fournaise. Within this context, we describe the implementation of the system, including the lava flow numerical modeling, the use of satellite sensor to derive the time-averaged discharge rate (TADR) and flow outlines, and the production and communication of lava hazard maps giving the probability of lava flow coverage during an on-going crisis. An initial map is shared within the first hours of any eruption and, when needed, the map is updated during the eruption. As we explain, these commonly understandable maps have been developed through open and iterative discussion between scientists, the observatory responsible for monitoring and reporting du- * https://reliefweb.int/report/democratic-republic-congo/ r-publique-d-mocratique-du-congo-strat-gie-de-la-r-ponseruption † Kīlauea Eruption Recovery. How much damage was caused by the 2018 eruption? https://recovery.hawaiicounty.gov/connect/faqs ‡ Copernicus Emergency Service mapping: https://emergency. copernicus.eu/mapping/list-of-components/EMSR546 § Cabildo de la Palma: https://riesgovolcanico-lapalma.hub.arcgis. com ties (Observatoire Volcanologique du Piton de la Fournaise Institut Physique du Globe de Paris; OVPF-IPGP), and local civil protection (État-Major de Zone et de Protection Civile de l'Océan Indien; EMZPCOI). We provide details on items that was requested by the EMZPCOI to improve the delivered map. We discuss the limitation of our protocol by providing details on the uncertainties of the model itself, the influence of the topography that needs to be regularly updated, and the uncertainties related to the estimation of the discharge rates from satellite sensors. Finally, we examine the challenges of implementing such a multi-agency protocol, including the timeline of production and advances for future crises, and consider the potential for possible export of our protocol to other effusive centers.

Background on mitigating the risks associated with lava flows
Risk mitigation associated with lava flow consists of preparedness, evacuation, and sometimes attempts to divert or delay the advance of the flow [for review see Harris 2015]. As highlighted by Harris [2015]: "preparation and planning for losses remains the best way to reduce the severity, or moderate the impact (of lava flow ingress into a populated area) while also alleviating the loss incurred".
During effusive eruptions, monitoring lava flow emplacement is a common action carried out by volcano observatories to track and assess land and infrastructure lost or at risk. Since the 1990-2000s, the development of event onset alerting tools using high-spatial resolution satellite sensors operating in the infrared have been providing support to observatories [Harris 2013;Harris et al. 2016]. More recently, numerical lava flow models (stochastic or deterministic) have also been developed to include satellite-derived source terms (mostly TADR in m 3 s −1 ) and to allow real-time lava flow emplacement simulations and projections. This has permitted construction of operational tools, such as generation of hazard maps to assess potential lava flow inundation area and to project lava flow paths [Wright et al. 2008;Ganci et al. 2012;Harris et al. 2019]. These tools have been provided as support for stakeholders involved in preparing and planning for loss (e.g. volcano observatory, civil protection, and/or national, regional, and local government agencies). Since 2007, the Hawaiian Volcano Observatory (HVO) has been evaluating the short-term threat to communities posed by lava flows through frequent airborne and satellite mapping of lava flow fields and assessing steepest descent paths [Kauahikaua 2007]. These techniques have regularly been used during effusive activity at Kīlauea to build hazard maps including likely lava flow paths, and routes for communication with emergency managers and the public have been welldeveloped Brantley et al. 2019]. In 2014, during the lava flow crisis in the vicinity of the town of Pāhoa, scientists from the HVO measured lava flow advance rates and provided a range of potential arrival times that were used to build emergency plans [Poland 2016]. During the 2018 VOLC V NIC V 5(2): 313-334. https://doi.org/10.30909/vol.05.02.313334 eruption in the Puna district, HVO rapidly produced preliminary lava flow path forecasts using the DOWNFLOW model of Favalli [2005], a commonly used stochastic model that computes the line of steepest descent (LoSD) and probable lava flow inundation area. Throughout the eruption, lava flow path simulations were run on regularly updated topography, and from active flow fronts, new fissures, and channel overflow locations, to yield maps of likely future flow directions and inundation areas [Neal et al. 2019]. These maps facilitated communication of up-to-date hazard information to emergency managers and provided situation awareness for eruption response field crews. On Mount Etna (Italy), the Istituto Nazionale di Geofisica e Vulcanologia (INGV) at the Etna Observatory (EO) models and tracks lava flow emplacement using a combination of satellite-derived discharge rate estimates and numerical models [Vicari et al. 2011]. INGV-EO employs the HOTSAT satellite monitoring system to convert infrared data from both MODIS and SEVIRI sensors into TADR that is then input in the numerical model GPUFLOW. GPUFLOW, an evolution of MAGFLOW, is based on a cellular automaton, physics-based model which can reproduce, spatially and temporally, the likely lava flow propagation [Vicari et al. 2011;Ganci et al. 2012;Del Negro et al. 2013;Bilotta et al. 2016]. The combination of HOTSAT and GPUFLOW is integrated into a web-GIS system, LAV@HAZARD, that produces real-time hazard assessment by providing time of propagation of lava flow fronts, maximum runout distance, and area of inundation ]. This operational tool was described and validated using a retrospective analysis of Etna's 2008-2009 flank eruption by Ganci et al. [2012]. It has since been operational by providing automatic reports including the radiant heat flux, TADRs, and a short-term hazard map. These products are communicated to INGV and local civil protection through periodic bulletins. Although the monitoring strategy implemented in LAV@HAZARD was designed for Etna, it has been adapted and employed during volcanic crises at other sites, such as Fogo at Capo Verde  and for Nabro in Eritrea [Ganci et al. 2020].
At Piton de la Fournaise, a well-established protocol between the OVPF-IPGP, local civil protection (EMZPCOI), and scientists from a multi-national array of institutes has allowed effective tracking of effusive crises and hazard management . Forecasting and modeling lava flow inundation areas through a satellite data-driven protocol to mitigate their impact started to be developed in 2014 ]. The DOWNFLOW model was used to provide the LoSD and probable lava flow inundation area [Favalli 2005] and FLOWGO was used to calculate the maximum runout of lava for a given effusion rate [Harris and Rowland 2001;Harris 2015]. This protocol and near real-time application, as tested and validated during the April 2018 eruption of Piton de la Fournaise, was presented in Harris et al. [2019] and its improvement and final, operational implementation is provided and discussed herein.

Eruptive activity at Piton de la Fournaise
Since the first documented eruption on La Reunion in 1708, 95 % of the eruptions of Piton de la Fournaise have oc-curred inside the Enclos Fouqué caldera ( Figure 1B) [Villeneuve and Bachèlery 2006;Chevrel et al. 2021]. The Enclos Fouqué caldera is open towards the ocean to the east and its formation has been accompanied by successive landslides [Bachèlery 1981], forming a horseshoe-shaped depression, hereafter called the Enclos. Vents generally open within the Enclos, with lava flows being contained within the depression and flowing downhill to the east and towards the ocean ( Figure 1B, C). Eruptions are fed by magma dikes or sills that propagate along preferential paths within, tangentially, or radially around the summit crater (i.e. the Dolomieu crater) of the central terminal shield inside the caldera and along three main rift zones orientated southeast, northeast, and 120°N [Bachèlery 1981;Dumont et al. 2022]. Since at least the last century, the majority of eruptions have consisted of fissurefed lava fountaining and lava flows .
Eruptive activity is controlled by cycles related to the evolution of the superficial edifice stress field that usually last from one to eleven years and control the location of eruptive fissures Got et al. 2013;Derrien 2019] and to magma source recharge [Vlastélic et al. 2018]. These cycles usually start with summit or summit-proximal eruptions and end with a distal or eccentric eruption that occurs more than 4 km from the summit and can be located inside or outside of the Enclos. In contrast, deep source-related cycles are one to three decades long and control eruption output rates and erupted volumes. These cycles end with summit crater collapse and voluminous eruptions (a "hotspot surge") followed by few years of quiescence [Vlastélic et al. 2018]. This was the case, for example, in April 2007 when the most voluminous low-elevation flank eruption so far observed on La Réunion ] marked the end of the 1998-2007 cycle [Vlastélic et al. 2018]. The April 2007 vent was located at 590 m a.s.l., and lasted one month (30 March-1 May, 2007) and produced 140 to 240 × 10 6 m 3 of lava Roult et al. 2012]. Lava flows were fed by a high effusion rate (around 200 m 3 s −1 ) and reached the island belt road, located 2.2 km from the vent, in less than 5 hours . By the end of the eruption, the belt road that links the south and north of the island had been buried over a length of 1.4 km by a 50 m thickness of lava ]. The April 2007 event was followed by three years of low intensity activity inside or close to the summit at a mean output rate of 1.8 × 10 6 m 3 yr −1 . This compared with 24 × 10 6 m 3 yr −1 towards the end of the 1998-2007 cycle [Roult et al. 2012] and was followed by four years of inactivity from 2010-2014 .
In June 2014, Piton de La Fournaise entered a new period of high eruptive frequency. The single day eruption on 20 June, 2014 was followed by four eruptions in 2015, four during 2016-2017, and then 14 during 2018-2021 ( Figure 1C). The 2014-2021 eruptive period was characterized by both short and long duration eruptions (0.7 to 55 days) and mean output rates of 1.5-18 m 3 s −1 (mean of 6.8 m 3 s −1 ; annual average of 14 × 10 6 m 3 yr −1 * ). The majority of eruptive fissures were located within 2.5 km of the summit above 1800 m a.s.l. and inside the Enclos ( Figure 1C). Exceptions were the two last eruptions of 2019 (starting on 11 August and 25 October, respectively) during which fissures opened in the Grandes Pentes (steep slopes) area at distances of 4.5 and 6 km from the summit at 1450 and 1000 m a.s.l., respectively. For both eruptions, lava flows reached the area of the Grand Brûlé (deeply burnt) and stopped 2 km and 220 m from the belt road during the two eruptions, respectively ( Figure 1C). Additionally, lava flows associated with the September 2016, July 2018, and February-March 2019 eruptions buried several hundred meters of the hiking trail network ].

Volcanic crisis management at Piton de la Fournaise
Although most eruptions take place inside the uninhabited Enclos, this area remains highly frequented, with at least a hundred thousand hikers access the summit every year [Derrien et al. 2018]. Furthermore, the lower eastern flank of the Enclos is crossed by the island belt road (RN2), which represents the only access between the southern to northern parts of the island along the east coast. Since 1950 the road has been cut 10 times by lava flows ]. On 12 occasions since the beginning of the 18 th century eruptions have fed lava flows outside of the Enclos (called Hors Enclos eruptions) where populations are at risk; this was the case recently, in 1977in , 1986. While the 1998 event produced short lava flows at a high elevation, the 1977 and 1986 eruptions fed lava flows that entered the populated coastal area. The April 1977 eruption required evacuation of Piton Sainte-Rose, burnt down or damaged~30 buildings (houses, the church, the police station and a gas station), and buried a portion of the village [Kieffer et al. 1977].
Following the 1977 eruption, the French Ministère des Universités, the département de La Réunion, the Institut National d'Astronomie et de Géophysique (now called CNRS) and the Institut de Physique du Globe de Paris (IPGP) established the Observatoire Volcanologique du Piton de la Fournaise (OVPF) with a remit to monitor the volcano. Additionally, a national response organization plan was established in the case of eruptions at Piton de la Fournaise by the French government for civil protection and associated agencies. This plan, still valid today, is named the "ORSEC-DSO (Organisation de la Réponse de SEcurité Civile-Dispositions Spécifiques Opérationnel), Volcan du Piton de la Fournaise", and defines the warning levels to be set, as well as the actions to be taken and restrictions to be implemented at each level, the protocol to change the level, and the communication strategy for alert diffusion . As part of the plan, in the case of a change in activity (i.e. onset of unrest, new fissure opening, end of eruption), OVPF must inform the civil protection headquarters of the Indian Ocean zone (EMZPCOI) and propose a change in the alert level. The decision to officially change the alert level is the responsibility of the Préfecture, (i.e. local government headquarters). The Préfecture then communicates with other actors (town councils, gendarmerie (police), central authorities, institutions, and the media).
The four levels of alerts are: • Vigilance: possible eruption in the medium term (a few days or weeks) or presence of risk (rockfalls, increased gas emissions, inactive but still hot and cooling lava flows ...); • Alert 1: probable or imminent eruption; • Alert 2: ongoing eruption; • Alert 2-1: ongoing eruption inside the Enclos without threat to the safety of people, property or the environment; • Alert 2-2: ongoing eruption inside the Enclos with threat to the safety of people, property or the environment; • Alert 2-3: ongoing eruption outside the Enclos with threat to the safety of people, property or the environment; • Sauvegarde: end of eruption.
As soon as Alert 1 is triggered, an inter-services meeting is organized at the Préfecture to (i) establish a situation report to inform the response services, and (ii) prepare the crisis management team and response system. The inter-service consortium of responders includes: • the police force (gendarmerie), • the town councils, • the fire department (SDIS: Service départemental d'incendie et de secours), • the department of environment, planning and housing (DEAL: Direction de l'environnement, de l'aménagement et du logement from Préfecture), • the national forest office (ONF: Office National des Forêts), and • if required, the departmental and regional council (conseil départemental de la Réunion, conseil Régional), the army (FAZSOI: Forces armées dans la zone sud de l'Océan Indien), the emergency health services (SAMU: Service d'Aide Médicale Urgente, and ARS: Agence Régionale Santé) are also included in the consortium.
From alert level 1 onwards, access to the Enclos Fouqué is strictly limited to authorized personnel.
During an eruptive crisis OVPF communicates daily with EMZPCOI via daily bulletins and phone calls. In the case of a "special situation", for example when there is an imminent threat to population or infrastructure, as is the case for eruptions located outside the Enclos or when there is a risk of the belt road being cut, the Centre des Operation de la Préfecture (COP) is activated. Activation of the COP means that the system to organize and implement mitigation actions for the population and infrastructure at threat, or road closure, is triggered. If the situation requires, one or more operational command posts (PCO: Poste de Commandement Opérationnel) and advanced medical stations (PMA: Poste Médical avancé), may also be installed as close as possible to the impacted area. This is a decision made by the Préfet (head of the Préfecture). To date, the COP has been activated only once, in April 2007.  At the alert level Sauvegarde, the eruption is over but restrictions on access to the volcano remain in place until activity indicators (seismicity, gas flux, and ground deformation) are considered to have returned to "normal". Reconnaissance is also carried out to check that the area is safe for reopening. The alert level then reverts to "vigilance".

LAVA FLOW NUMERICAL MODELING AT OVPF
A number of numerical models exist to support lava flow modeling (see review by Dietterich et al. [2017]). Here, we combine the stochastic model DOWNFLOW [Favalli 2005] with the deterministic model FLOWGO [Harris and Rowland 2001;Harris 2015] to support hazard assessment at Piton de la Fournaise.
DOWNFLOW provides the most likely lava flow paths, including the LoSD and area of coverage. For responding to current crises at Piton de la Fournaise, this model was calibrated by fitting the output flow coverage to the actual areas of all flow fields since 2016 .
FLOWGO calculates the runout distance of a lava flow along a slope line for a given effusion rate [Harris and Rowland 2001]. It is 1-D model adapted for cooling-limited lava flow in a channel of uniform depth that, once calibrated with a suitable rheological model, only needs the slope from the vent along the steepest descent line and a discharge rate as its source terms (see Chevrel et al. [2018] for more details). Once initialized, to assess flow runout, FLOWGO takes a control volume of lava and tracks its cooling and crystallization down-channel until it becomes too viscous for further forward motion to occur, i.e. velocity approaches zero [Harris and Rowland 2001]. FLOWGO was calibrated to be run on the LoSD obtained via DOWNFLOW and using a retrospective analysis of the April 2018 eruptions ]. If needed, new calibration can be carried out at any time if source terms are known to change or if the model fit with the actual flow begins to break down. It must be stressed that FLOWGO output is only appropriate for channel-fed, cooling-limited flow of basaltic composition.

Original chain of tasks
The rapid response lava flow simulation protocol began with the eruption of June 2014 ]. Between 2014 and 2018, the chain of tasks comprising this protocol involved multinational partners and each step was operated by a different institute in a different country and was collated on a centralized reporting page ]. As part of this protocol, five sequential actions were executed: 1) Action 1: At the start of an eruption, OVPF communicated the vent or fissure coordinates (obtained either via visible observation or location of the tremor that is available online on the WebObs platform ) to INGV in Pisa (Italy) where a lava flow inundation map was produced using DOWNFLOW on the most recent Digital Elevation Model (DEM).
2) Action 2: The DOWNFLOW-derived LoSD was sent to the Laboratoire Magmas et Volcans (LMV) at the Université Clermont-Auvergne (UCA, France) where the FLOWGO model was run to estimate maximum the lava flow runout distances.
Note that uncertainties for runout distances are of the order of 30 % .
3) Action 3: The key source term for FLOWGO is the discharge rate. This value was derived from satellite-based thermal infrared images with the MIROVA platform at the Universita di Torino in Italy * and with the HOTVOLC platform at the Observatoire de physique du globe de Clermont-Ferrand (OPGC) at UCA † two volcanic hotspot detection systems. Uncertainties on the satellite-derived TADRs are of the order of 50 % for HOTVOLC and 35 % for MIROVA (see discussion). The cumulative volume can also be calculated by integrating TADR values over the period of time in between two satellite images using the trapezium rule. 4) Action 4: The resulting lava flow area coverage obtained via Action 1, and the lava flow runout maximum distance obtained for the given discharge rate (Actions 2 and 3) represented as a graph of lava velocity versus distance from the vent to the flow front were then communicated to OVPF ]. 5) Action 5: In parallel, the evolution of lava flow field (vent location, channel dimension, flow outline) was obtained from airborne photogrammetry acquired by OVPF. In addition, interferograms and coherence maps derived from InSAR via the OI 2 platform developed at the OPGC, UCA ‡ and satellite infrared and visible imagery acquired by the ASTER Urgent Request Protocol [Ramsey et al. 2016] were used to assess zones covered by lava. These products were also used to close the loop by up-dating and cross-checking the performance of DOWNFLOW and FLOWGO ].
This protocol was triggered in real time at the start of any given eruption through release of an email from the director of OVPF with vent coordinates and was designed so that it could be repeated and updated as the eruption conditions evolved (e.g. new vents opening, extension of tubes, change in TADR). The response protocol was reviewed and updated throughout each of the ten eruptions from June 2014 to April 2018, at which point the protocol became fully operational ].

Improved chain of tasks
The system became centralized in 2019 to improve efficiency. The execution of DOWNFLOW and FLOWGO was combined in DOWNFLOWGO as part of the upgrade Peltier et al. 2021], thereby merging Actions 1 and 2. This was possible through PyFLOWGO, an integrated python code that allows the actions to be executed by a single operator at a single location on any operating system . Action 3 was also streamlined because the original process required a new simulation to generate the flow runout projection with each change in the TADR, which was time-consuming and inefficient. Simulations are now run for a range of likely effusion rates (typically between 5 and 50 m 3 s −1 at increments of * http://www.mirovaweb.it/; [Coppola et al. 2016]. † https://hotvolc.opgc.fr; [Gouhier et al. 2016]. ‡ https://opgc.uca.fr/volcanologie/oi2; [Hrysiewicz 2019;Richter and Froger 2020]. 5 m 3 s −1 ). In essence, once the vent coordinates are entered in the code line, a raster of the most probable area to be covered, the path of steepest descent, and the runout locations (coordinates and elevation) along this path are computed. This is carried out on a 10 m resolution DEM (here we use the DEM acquired via LiDAR in 2010) that is 130 × 10 6 m 2 and where 10,000 DOWNFLOW runs are calculated, as calibrated by Chevrel et al. [2021]. Then runs of FLOWGO for effusion rates, at 10 m steps down the LoSD are made, this being the ideal iteration step size to ensure numerical convergence ] to give the runout points. The complete computing time is less than two minutes. The results are then imported into Q-GIS, a free and open-source Geographic Information System, that already contains a pre-prepared template having all the needed layers, including locations of the OVPF monitoring network and all features required by civil protection (see Section 5 for details). This improved protocol is given in Figure 2 and only requires the vent coordinates of the starting eruption to prepare and deliver the map, with the map being prepared in a few minutes by a single operator.
The delivered short-term hazard map (e.g. Figures 3, 4, 5, 6 and 7) provides: 1) The path of steepest descent from the main vent location (in case of multiple vents along a fissure, several simulations may be run); 2) The probable area of lava flow coverage (i.e. the probability that a pixel will be hit by lava); and 3) The runout distances (i.e. where the lava should stop) for the pre-defined effusion rate range along the LoSD.
The map is then used as a guide to assess probable runout given current TADR obtained from the satellite detection systems. Actions 3 and 5 are still implemented through email exchange and/or consultation the dedicated websites. This includes update of TADR from MIROVA or HOTVOLC to allow the look-up-based assessment to be revised, as well as provision of ASTER images and InSAR coherence maps to allow on-the-fly validation of the hazard map and possible update.
The map is first reviewed and approved by the OVPF director, who is then responsible for transmitting the map to the local civil protection. Due to difficulties of representation, uncertainties are not represented on the maps but known by all actors. For this, regular explanation of the modeling is done during meetings with civil protection (see Sections 5 and 6). This is essential to ensure that all duty staff know that these maps are a first-order estimation of the area to be covered and that the lava can always take another direction or go further than predicted. By common agreement, the maps are not open to the general public during an eruption. This is to avoid misinterpretation from untrained operators and the spread of false information. However, following any given eruption, the maps are published in the OVPF monthly bulletin § . § (e.g. The third level (dark grey) are the two employed software DOWNFLOWGO and any Geographic Information System. DOWNFLOWGO is composed of (i) DOWNFLOW that needs the topography (DEM) and the calibrated simulation parameters (Δh and N) and computes the line of steepest descent (LoSD) and the probable lava flow area, and (ii) PyFLOWGO that uses the lava thermo-rheological properties to compute the runout distances along the LoSD for a given effusion rate. The GIS needs to contain all required layers (geology and infrastructure) and a map template so to produce the consistent final product that is the short-term lava flow hazard map to be delivered to the authorities. Acquired lava flow outline by field and airborne survey or satellite data acquired during the eruption can be added as a layer as and when available.

STUDY CASES: SELECTED ERUPTIONS
To illustrate our operational near real-time protocol, we present the response to three eruptive events as case studies, these being the eruptions of 25-27 October, 2019 (an eruption that threatened the road belt), 10-16 February, 2020 (an eruption that started only 20 minutes after the start of the seismic crisis), and 7-8 December, 2020 (when the response protocol was triggered in a record time after the eruption onset). Additionally, we provide a hypothetical case study of an eruption outside of the main caldera in an inhabited area, based on the scenario proposed in Tadini et al. [2022]. For each case study we describe the eruption timeline (given in the local time zone, i.e. RET, UTC+4) and include the maps that were communicated to the authorities (note that original maps have legends in French but they are here translated into English), a map showing the final lava flow outline, TADRs as function of time as obtained from HOTVOLC and MIROVA, and the respective cumulative volume. The end of the eruption, as declared by OVPF when tremor ceases, is also given.

25-27 October, 2019
The 25-27 October, 2019 eruption was characterized by the location of the vent that was the lowest in elevation since the 2007 eruption, and the island belt road was threaten by the lava flow. A seismic crisis started at 04:15 RET (UTC+4) on 25 October, 2019 after~14 days of edifice inflation. At Piton de la Fournaise, such seismic crises have been interpreted as a sign that the shallow magma reservoir has pressurized to the extent that a diking event occurs, triggering magma propagation out of the chamber and towards the surface . After less than 1 hour of seismic crisis, the location of the earthquakes and of the origin of ground deformation indicated magma propagation towards the eastern flank. A total of 827 shallow volcano-tectonic earthquakes below the summit and ground deformation of~12 cm were recorded over the following ten hours [OVPF 2019]. Volcanic tremor, which at Piton de la Fournaise is associated with magma reaching the surface with opening of a fissure, and thus interpreted as the onset of eruption, appeared on the seismic records at 14:40. From the location of the tremor source, the eruptive site was located at a low elevation in the Southeast sector of the Enclos in the Grandes Pentes area ( Figure 1). We produced the first hazard map at 16:00 local time using this approximate location ( Figure 3A). This initial product was communicated with OVPF only because the exact location of the eruptive fissure was not yet confirmed through visual observation, and was used to assess what parts of the monitoring network could be at risk and what likely scenarios to prepare operations for.
An aerial survey via helicopter took place 2.5 hours after the eruption onset and revealed two eruptive fissures in the southsoutheastern sector of the Enclos, at 1060 m and 990 m a.s.l., respectively. At the time of the survey, a third, smaller upper eruptive fissure was no longer active, and activity had focused on the lower fissure. At 16:45, a first TADR of 6 m 3 s −1 was obtained from HOTVOLC. This value was likely underestimated due to extremely cloudy conditions and was thus treated with caution. From the improved vent location obtained from observations on the overflight, the lava flow simulation was rerun and a new map produced ( Figure 3B), which was sent to EMPZCOI at around 18:00 (3.3 hours after the start of the eruption). From this map it became clear that the belt road could be cut given lava flow fed at any discharge rate greater than 20 m 3 s −1 .
At 22:30, lava discharge rate from HOTVOLC exceeded 25 m 3 s −1 but a first MODIS image was acquired 3 hours later at 01:30 (on 26/10, 11 hours after the start of the eruption) from which a lower discharge rate of between 7 and 14 m 3 s −1 was obtained. By the next morning, HOTVOLC indicated a decline in discharge rate with values <10 m 3 s −1 (Figure 3D), thus suggesting also a decrease of runout potential. On the same day (26 October), improved weather conditions also allowed a field survey to obtain a more accurate vent location. The hazard map was subsequently updated and sent to EMPZCOI at the end of the day. This map showed the exact location of the main vent and the position of the lava flow front at that time, which matched the simulated lava flow path and runout for a TADR of 20 m 3 s −1 ( Figure 3C).
The eruption stopped the following day, on 27 October at 16:30 RET after 1 hour of gas pistoning activity [OVPF 2019]. Overall, lava discharge rates estimated from the satellite platforms were up to 27 m 3 s −1 for the first 20 hours of the eruption, which then decreased to <10 m 3 s −1 thereafter (Figure 4A), although these values were likely underestimated due to extensive cloud cover. Final cumulative volume obtained from TADRs was 1.0±0.35 × 10 6 m 3 and 0.8±0.4 × 10 6 m 3 for MIROVA and HOTVOLC, respectively ( Figure 4B), which represents 50 to 62 % less than the total lava volume of 1.6 × 10 6 m 3 obtained from photogrammetry (for details on the method, see Derrien [2019]), which confirms an underestimation in TADR.
The final lava flow length was 3.6 km and the flow front stopped 230 m short of the belt road. The flow front location corresponded to the flow front runout modeled for a TADR of 20 m 3 s −1 ( Figure 3D). The final lava flow field outline from post-eruption airborne photogrammetry revealed a good fit with the simulated path ( Figure 3C), further showing that our calibrated model performed well. Note that uncertainties on runout distances are of the order of 30 % , which in this case is ±1 km. Given this degree of uncertainty it was thus not possible to disregard the threat to the road.

10-16 February, 2020
This eruption was marked by a very short time between the onset of precursory activity and the start of the eruption. Indeed, the eruption began only 20 minutes after the onset of  the seismic crises which gave little reaction time to produce the short-term hazard map.
On 10 February, 2020, a seismic crisis started at 10:27 RET (06:27 UTC) with more than 200 shallow volcano-tectonic earthquakes in 30 minutes. This was accompanied by rapid surface inflation of up to 27 cm [OVPF 2020]. Volcanic tremor was detected at 10:50, located on the southeastern flank of the summit shield. Due to fog, no information on the precise fissure location was known until 11:15, when a local resident called the OVPF while observing a gas plume sourced on the east of the volcano, in the vicinity of an existing scoria cone (Marco Crater). The first simulation was therefore run arbitrarily from a few hundreds of meters above of this cone and a map was delivered to EMPZCOI at 12:00 ( Figure 5A). Given that the fissure extended across the boundary of two municipalities, EMPZCOI requested that these boundaries be added to the map. The first lava discharge rate estimation with MIROVA was obtained at 13:20 (2.5 hours after the onset of the eruption) and indicated a range of 3 to 5 m 3 s −1 . HOTVOLC provided an initial, isolated (single-point) peak at 36 m 3 s −1 , after which values were more consistent with those of MIROVA. Due to the cloudy conditions these values were likely underestimated, but these initial estimates allowed the first assessments of likely flow runout ( Figure 5B).
In the early afternoon, an aerial survey revealed that the main eruptive fissure was oriented parallel to the slope from west to east, extending downslope between 2450 m to 1990 m a.s.l. It was not possible to place a single point for the source of emission because the fissure was parallel to the slope. We therefore performed a series of simulations along the fissure to consider the full range of source elevations. To avoid overcrowding the map, we decided to only represent the paths from the upper and lower ends of the fissure, these representing the en-member bounds on flow runout. An updated map was then delivered to EMPZCOI at 15:00 on which MIROVA-derived TADRs were highlighted ( Figure 5B). For this new map, EMPZCOI asked for the addition of vegetated areas so that they could quickly inform the relevant municipality as where their respective fire departments should intervene in case of wildfire. Our simulations suggested that the lava flow would reach the Grande Pentes and the Grand Brulé at discharge rates of about 20 and 30 m 3 s −1 , respectively, and that the municipality of Sainte Rose, to the north, was the most likely to be affected by wildfires ( Figure 5B). Two other smaller flows located to the southwest of the main fissure were also observed but these were not active by the time of the overflight.
On 11 February, a coherence map was acquired via InSAR and the lava flow field outline was extracted by the OI 2 platform. This outline showed a lava flow front near the expected runout distances for discharge rates of 10-12 m 3 s −1 , meaning that the flow would have already reached its cooling-limited length if the discharge rate was as reported by HOTVOLC and MIROVA ( Figure 5C). For the next two days, discharge rate estimation was not possible due to cloud conditions. Images were acquired on February 13 and 14 (63 to 109 h after the start of the eruption), indicating TADRs of 8-13 m 3 s −1 , as estimated by MIROVA and around 5 m 3 s −1 as estimated by HOTVOLC (Figure 5D).
Seismic tremor decreased abruptly around 14:00 on 15 February, and the eruption stopped 24 hours later at 14:12 on 16 February. Once the final lava flow outline had been obtained from photogrammetry after the eruption, a new simulation was performed from the exact location of the main lava source to compare performance of the modeling ( Figure 5C). This revealed that the flow had indeed reached its coolinglimited length within the first two days, and did not extend further; with the eruption continuing for three more days at approximately the same discharge rate. During this period, instead of extending downslope, the flow extended laterally, forming branches to the north as revealed by the difference between the outline obtained on 11 February by InSAR and final flow outline ( Figure 5C). Thickening of the flow field by superposition of lava units could have also taken place. In total, the erupted volume was estimated from the aerial photogrammetry at 3.8 × 10 6 m 3 , and this was consistent with the value of 3.5±1.2 × 10 6 m 3 obtained with MIROVA ( Figure 5D). The HOTVOLC system provided a lower volume of 2.2 ±1.1 × 10 6 m 3 , and we note that this underestimation was likely due to cloudy conditions.

7-8 December, 2020
For the 7-8 December, 2020 eruption the response protocol was triggered in a record time, with the map being provided to EMPZCOI just a few minutes after the initial aerial survey of the eruption site, which itself was just a couple of hours after the eruption began. On 4 December, 2020, after a slight ramp-up in seismicity starting on 3 December and inflation in the first few days of December, a dramatic increase in seismicity began at 05:10 RET and was accompanied by small levels of ground deformation (<1 cm) below the Dolomieu crater. Although this crisis did not lead to an eruption right away, the seismicity remained high for three days until a second crisis began on 7 December at 02:28. This new crisis was accompanied by larger ground displacement (reaching about 30 cm) and localized to the south-southwestern flank of the terminal shield. Two hours later, at around 04:40, the volcanic tremor was recorded and localized on the southwestern sector of the Enclos [OVPF 2021].
A first overflight of the eruptive site was possible at 07:30, which allowed OVPF to precisely locate the eruptive fissures just after the sunrise. On return to the observatory and in less than ten minutes, the coordinates of the eruptive site were entered and the first lava flow hazard map was produced, approved, and transmitted to EMZPCOI, by 08:00 ( Figure 6A). The production and delivery of the map was extremely rapid because all operators were onsite at OVPF and because the map template had been pre-prepared by setting up all the GIS layers over the previous days knowing an eruption was imminent. The eruption was located near the heavily vegetated south wall of the Enclos, which also necessitated rapid notification to EMZPCOI due to the possibility of wildfires. Discharge rate estimates given by HOTVOLC during the first hours of the eruption ranged between 10 and 42 m 3 s −1 , and, given this large range, it was not possible to discount lava reaching the vegetated rampart. The first map communicated to EMZPCOI therefore did not contain the lava flow runout estimates, but only the lava flow path and location of the front position at the time of the flight survey ( Figure 6A). At 10:05 on the same day (7 December), MIROVA revealed a TADR of 18±6m 3 s −1 . These agreed with the HOTVOLC data for the same time ( Figure 6B). A new map was thus produced and communicated to EMZPCOI showing that, with this lower value, it was unlikely that the lava flow would reach the caldera wall ( Figure 6C). However, on the map it was clear where the point of contact would occur if the lava discharge rate increased. By the end of 7 December, discharge rate had decreased to <2 m 3 s −1 and was accompanied by a gradual decrease in tremor. The eruption stopped the following day (8 December) at around 07:15 after 26 hours and 35 minutes of effusive activity. The post-eruption lava flow field outline from aerial photogrammetry showed that flows extended to 2 km, with the flow front stopping 500 m short of the caldera wall. A Sentinel-2 image (acquired at 10:35 RET on 7 December and processed the next day) also provided the lava flow field area. The fit between the flow field outline and the model was good in the proximal zone but diverged in the distal part ( Figure 5D). This discrepancy was due to the presence of a previous lava flow (of April 2018) that was not present on the DEM. As discussed later, the DEM needs to be regularly updated at such an active site, although unfortunately an update could not have been done in time for this eruption. The total erupted lava volume was estimated at 0.8 × 10 6 m 3 from photogrammetry, which was consistent with the cumulative volumes obtained from discharge rate time series (0.8±0.3 × 10 6 m 3 and 0.9±0.45 × 10 6 m 3 , with MIROVA and HOTVOLC respectively).

Hypothetical eruption outside of the caldera
This last study case focuses on a hypothetical eruption that would occur outside of the main caldera (i.e. Hors Enclos eruption) in an inhabited area. This scenario is based on the expert elicitation exercise proposed by Tadini et al. [2022]. We chose this case to test our protocol for an eruption that could potentially impact population in inhabited areas.
In this hypothetical scenario, tremor started at 12:05 RET on Day 0 on the north-western flank of the volcano. Eruptive fissures opened about 15 km north-west of the caldera in the town of La Plaine des Palmistes. Initial thermal anomalies were detected at the same time, with an estimated TADR of 120 m 3 s −1 . Such a high value is consistent with the actual highest discharge rate to date of 129 m 3 s −1 recorded by HOTVOLC on 24 August, 2015, at 15:45 on Piton de la Fournaise.
A vent opening in a populated zone requires particularly rapid response, and our previous case study in Section 4.3 (7-8 December, 2020 eruption) shows that, if all actors are together on site, such map can be produced and sent to EMZPCOI at best within 10 minutes of the communication of the vent location. In the case of the expert elicitation exercise proposed by [Tadini et al. 2022] only the path of steepest descent was provided to the experts and this was about 6 hours after the start of the eruption because it included reports from initial field reconnaissance by OVPF-IPGP staff, confirming vent locations and flow front locations. Here we show the first hazard map that we could have been produced following our protocol ( Figure 7A). From this map it is possible to rapidly identify the potentially affected infrastructure. Note that simulations do not capture interaction with buildings that could possibly slow down or locally divert the lava flow [e.g. Dietterich et al. 2015]. Proximally, we see the flow impacting the north side of La Plaine des Palmistes ( Figure 7B) and distally the predicted range of potential lava flow paths is confined to a deep ravine, except where the flow enters the ocean. Such ravines are characteristic of the geomorphology of La Réunion, and in the case of an Hors Enclos eruption, could be advantageous because these features may confine the lava. Note, however, that in some areas of the Piton de la Fournaise the Digital Terrain Model (DTM) is not reliable due to dense tropical vegetation. The use of a DEM calculated from the DTM may therefore lack vertical precision because the canopy tends to smooth out the details of the terrain [Stevens et al. 1999;Bigdeli et al. 2018]. The importance of the gullies in channelizing the lava flow should therefore be treated with care when lava flow simulations are performed in such dense tropical forest areas. Such a situation happened in Hawai ' i during the Pāhoa lava flow of 2014, where the lava flowing through a dense forest entered a preexisting series of ground cracks [Poland 2016], so the progress of the lava flow could only be tracked following the gas plume .
Because the eruptive fissure is situated at high elevation and about 15 km from the ocean, the modeled runout distances indicate that at the discharge rate given in this scenario (120 m 3 s −1 ) the lava would not reach the ocean but would require even higher discharge rate (>250 m 3 s −1 ) in order to enter the ocean ( Figure 7A). In such a case, the flow would cut all towns and roads between La Plaine des Palmistes and the populated coastal zone around Saint-Benoit. Now, if we compare with an analogous large Piton de la Fournaise eruption of such high discharge rate, such as that of 2007, lava tubes could form [Rhéty et al. 2017]. Formation of a tube will greatly reduce cooling rates by enhancing flow insulation, allowing the flow system to extend much further at any given discharge rate [Keszthelyi and Self 1998]. Considering that a tube system would develop, we can simulate a tube condition by assigning the crust cover fraction in FLOWGO to be 100 % and a low temperature [Rowland et al. 2005]. Applying this, we find that the lava could travel from the vent to the ocean (at a distance of 15 km) at discharge rates as low as 35 m 3 s −1 ( Figure 7A).

IMPROVING THE DELIVERED HAZARD MAP
Between 2014 and 2018, the products shared with OVPF comprised a map with the DOWNFLOW results (showing the probability distribution of lava flow paths) and the results of FLOWGO (presented as a graph showing the lava velocity versus distance from the vent to the flow front for the range of MIROVAgiven TADR) ]. The OVPF director could then inform EMZPCOI of the possible runout of the lava flow, but no physical product was shared. The primary format of communication was a telephone exchange between the OVPF and EMZPCOI, supported by information from visual observations made during flight surveys and/or field crews. Because the communication lacked product support, officers in charge of civil protection actions at the Préfecture took much  Palmistes. This map shows the runout distance for a given discharge rate in the case of an open channel system (as usually first assumed, in blue) and in the case formation of a tube system (in green). [B] Zoom of the map to identify proximal buildings and roads that will be the first to be at risk. caution by, for example, over-estimating the zone at risk in their assessments areas to be closed. The belt road to the east of the Enclos would also be closed when lava flows arrived in the Grandes Pentes area at a distance of 3 km from the road. The other difficulty with this communication format was, without a hazard map, assessment of the situation consisted of simply pointing to the potentially affected area on a topographic map without precision of the probability distribution of lava flow paths. A first short-term hazard map was built for the April 2018 eruption ] but the first time a short-term hazard map was shared with EMPZ-COI was for the August 2019 eruption. Since then, maps have been generated and shared for all seven eruptions between 2019 and 2021, with maps being sent to EMZPCOI within the first minutes to a few hours of eruption onset. These dedicated hazard maps have allowed EMZPCOI to have a better sense of the likely flow inundation area, front position and advance of the lava flow with discharge rate. Most importantly, the maps have allowed EMZPCOI to plan their response actions with higher degrees of confidence (e.g. to decide when and where to close the road, where to deploy wildfire countermeasures, or even when and where to evacuate populations, including hikers, and/or to close tourist overlooks that are identified as potential fire traps). Throughout the years the maps have been improved after each eruption upon specific requests from the EMPZCOI. The first short-term hazard maps only included lava flow path and runout distances tailored to the vent location and discharge rate of the specific eruption in hand and this was presented on a greyscale hill-shaded DEM background ( Figure 8A). Now, and as detailed successively, maps include a more complete background where remarkable items and vulnerable structures are visible ( Figure 8B). To ensure that the maps are as useful as possible in aiding EMZPCOI in improved crises management and that communication within EMZPCOI as well as with OVPF was fully informed, several exchanges and meetings were organized between EMZPCOI and OVPF, as well as LMV-OPGC. These exchanges served as a means of training for the emergency managers, so that EMZPCOI fully understood the scientific approach behind the modeling, including uncertainties in models. In return, it allowed the scientific group to understand the needs of the emergency managers. This work focused on improving the maps in terms of graphical representation following constructive comment from EMZPCOI. As a result, there is a great deal of difference between the first map produced for the April 2018 eruption and one of the most recent maps generated to date for the April 2021 eruption (Figure 8).
Map changes as a result of user needs include: • Color of lava inundation: One of the first exchanges was related to the colors. The blue shades, as presented in the first map ( Figure 8A) were too confusing as, for EMZCOI, the colors related to water flooding rather than lava inundation. Thus, as recommended by Thompson et al. [2015], colors "commonly associated with volcanoes and hazard, and [which] reduce the potential for confusion with hydrological hazard map types" were used (see also Harris et al. [2012]). This in-volved application of a yellow-to-red sequential color scheme, with increasing the probability of inundation ( Figure 8B).
• Boundaries of municipalities: Drawing the boundary of the two municipalities that shares the Enclos (Sainte-Rose to the north and Saint-Philippe to the south) was requested so that the right municipality could be quickly inform if needed. This is fundamental in cases where the lava flow reaches the road, as Sainte-Rose is responsible for closure, diversion and replacement/repair to the north, and Saint-Philippe to the south. It also helps the two municipalities to prepare for losses. To fulfil the request, the municipality boundaries where simply added to the background ( Figure 5A).
• Land cover: Including the vegetated areas was requested to visualize where potential wildfires could start. This allows civil protection to alert the local fire department(s) so they could intervene in time to stop fires igniting or spreading. This was also of use to the French National Forestry Office in case any endemic flora and fauna at risk could be saved. To complies with the request, we therefore added the vegetated areas as a GIS layer to the greyscale hill-shaded DEM ( Figure 5B).
• Well-known features: To improve easy geographical location of the potential lava flow impact, the map needed to include well-known features, such as trails, past lava flows, and well-known scoria cones. The anonymous shaded-relief background which is suitable for scientific purposes and model testing, turn out to not be appropriated for outreach and communication purposes. We therefore finally came into the conclusion that the most meaningful and usable background for EMZPCOI was the commonly used map of the sector from the Institut national de l'information geographique et forestiere (IGN) ( Figure 8B). This provides real added value to the product for response purposes and allows for rapid decision-making by the authorities.
• Specific items on the maps: EMZPCOI wanted the possibility to see specific items on the maps as organized into operational, strategic, and human resources. Operational resources include vital centers (rescue centers, hospitals and clinics, military camps, town halls, etc.) and strategic sites (airport, water reserves, gas stations, etc.). Strategic resources are roads, some specific sites and establishments such as schools and retirement homes, as well as electricity, phone, and internet networks. Human resources comprise homes, schools and colleges, as well as commercial or industrial areas. EMZP-COI therefore provided these layers as GIS-compatible format so that they could be implemented within our map template. Upon request they are activated to be visible or not.
The maps will certainly continue to be improved in the future. Indeed, some requests have been made but have not yet been fulfilled. For example, a critical tool for civil protection issues would be to automatically assess the number and location of people and buildings likely to be affected by the eruption, and therefore allow improved evacuation planning and loss preparation. This could potentially be done by adding the information of the map or providing a separate map where expected losses could be highlighted. Because (since the production of the maps) inhabited areas have not been threatened, this has not yet been implemented. Another request from the EMZPCOI is to include an assessment of the time before the lava flow front would reach critical morphological features (such as the caldera wall, the Grandes Pentes, or the coast) or infrastructure (such as the road and inhabited areas). This relates to the improvement of the the numerical model itself, rather than just the map; accordingly, this has not yet been provided.

Uncertainties of the modeling and future improvements
Between 2019 and 2021, ten eruption responses at Piton de la Fournaise used our lava flow response protocol. Because prior work Rhéty et al. 2017;Harris et al. 2019;Chevrel et al. 2021] had gone into initialization and validation of the base model (DOWNFLOWGO), our confidence in the numerical model was sufficient to ensure a fast response and evolution of the map product. However, to ensure a complete understanding of the delivered map, it is crucial to communicate the uncertainties and the limitations of the model to the end-user. Here we first present the uncertainties related to the calibration of the model then we discuss the uncertainties due to of the model itself and finally the potential model mismatch due to the topography.
For DOWNFLOW, a key uncertainty is that the accuracy of the projected lava flow area depends on the calibration of the numerical code that computes the steepest descent paths, while randomly applying a given vertical perturbation (Δℎ) at every pixel of the DEM during each of iterations [Favalli 2005]. The code computes a spectrum of possible paths of the flow from the vent to the edge of the DEM by iterating over a large number of runs ( ), producing a probabilistic estimation of the lava flow inundation area. Calibrating DOWNFLOW consists of finding the parameters Δℎ and that are able to best fit the model output to the actual lava flow area. For Piton de la Fournaise, Chevrel et al. [2021] found that these parameters may change with the volcano eruptive cycle and the available DEM resolution. For lava flows emitted between 1998 and 2007, simulations on a DEM produced in 1997 that has a 25 m resolution are the best for > 6,000 and Δℎ > 4 m. While for lava flows emplaced between 2010 and 2019 on a DEM from 2010 with a 10 m resolution, the best fit parameters are = 10,000 and Δℎ = 2 m. The best fit parameters changed because of the increasing DEM resolution and because after 2007, the lava flows are thinner and shorter ]. This last calibration ensures that more than 50 % of the real lava flow area is well-considered by the simulation ]. These best fit parameters are therefore used in the protocol, but must verified at each eruption to avoid over-or under-estimation of the lava flow coverage.
Changes in the initialization parameters for FLOWGO correspond to changes in lava composition, temperature, gas content, and crystallinity, all of which affect the rheology of the lava and thus also the flow velocity, thickness, and length. Harris et al. [2019] showed that variation in crystal and bubble content could affect the lava flow length by 28 %, while uncertainty on surface temperature could affect the lava flow length by 23 %. Moreover, Thompson and Ramsey [2021] showed that variable emissivity may have a measurable effect on heat flux (>30 % decrease), which translates to potentially of longer flows (~7 % increase). In addition,  argued that a two-component emissivity model produces flows that are a maximum of 34 % longer than flows modeled using a single constant emissivity. Here, we use the same input parameters for FLOWGO as defined by best-fitting to lava channel dimensional, crystal, and temperature properties for well-constrained channel-fed flow at Piton de la Fournaise by Harris [2015], Harris et al. [2017], Rhéty et al. [2017], and Har-VOLC V NIC V 5(2): 313-334. https://doi.org/10.30909/vol.05.02.313334 ris et al. [2019]. Based on post-event cross-checks of the model performance and the fact that lava composition has not been highly variable over the considered period, the model input parameters have been proven to still be valid ( Figure 3D, 5C, 6D). However, should the validity begin to break down we are able easily to change these parameters; moreover, the best fit is always appraised as part of model debrief at the end of each eruption. The modeled lava flow runout length uncertainty is 30 %. This is not yet represented on the map but this is regularly communicated to the end-user. Future improvement could include this uncertainty directly on the map.
Another source of modelling inaccuracy originates from the model itself. DOWNFLOWGO combines a range of probable lava flow paths and provides the cooling-limited distance the lava could reach down the LoSD. The model is unable to simulate levées, breakouts, overflows, lobes, tubes, bifurcation, etc. [Tarquini and Favalli 2011]. Also, FLOWGO is a 1-D model and assumes that the lava flow is channelized and cooling-limited [Harris and Rowland 2001;Rowland et al. 2005]. Collectively, this means that complex lava flow fields, volume-limited lava flows, and long duration flows cannot be modeled accurately because syn-eruptive topography changes that would influence the flow are not captured on the DEM. Nonetheless, at Piton de la Fournaise, effusive events are of relatively short duration. Additionally, this numerical model is computationally inexpensive, rapid, and can be operated on any operating system. Therefore, DOWNFLOWGO is, so far, the best option to rapidly build short-term hazard maps that show the most probable inundation area and the runout length for the given discharge rate at the beginning of an eruption when flow is in a well-defined channel. In the event of a long-duration flow, update of the DOWNFLOWGO modeling can be operated at the site of breakouts ] and insulation condition changes from open channel to tubes can also be modeled (see Section 4.4). Future options to better incorporate syn-eruptive lava flow morphology and flow front propagation could include use of a physics-based model, similar to MAGFLOW [Del Negro et al. 2013], which can spatially and temporally reproduce lava flow propagation [Cappello et al. 2022]. At Piton de la Fournaise, the Rheolef model has been used [Bernabeu et al. 2016], and VOLCFLOW [Kelfoun and Vargas 2016] has been tested on previous eruptions [Lemaire 2022]. Further work is needed to validate the use of such models during an eruption crisis and to ensure that they can be implemented in a timely fashion.
Errors in the modeled lava flow path may also come from error on the terrain topography [Tarquini and Favalli 2010]. In the Enclos, the high frequency of eruptive activity results in frequent change of the topography and these changes cannot always be added to the available DEM and therefore will not be considered by the model. For example, in this study we showed that the simulation of the lava flow emitted during the eruption of February 2020 did not compare well with the distal part of the flow. Here, the lava flow divided into two branches, while the modeled steepest descent line ran between the two branches ( Figure 6). This discrepancy is due to the April 2018 lava flow that was not present in the DEM that we used. This shows how important it is to update the DEM after every erup-tion. At Piton de la Fournaise, the frequent bad weather conditions on the volcano flanks and difficult terrain prevent the use of Unmanned Aerial Vehicles (UAV) for rapidly producing DEMs. Although this has been done in some cases and has been being completed by stereophotogrammetry from aircraft surveys [Derrien et al. 2015;Derrien 2019;Derrien et al. 2020], update of the DEM has not been systematic. DEMs are also produced by space agency services from optical sensors such as ASTER GDEM or radar sensors such as SRTM and tanDEM-X [Zink et al. 2014]. Additionally, at Piton de la Fournaise, with the help of the Centre National d'Etudes Spatiales (Kalidéos program), tri-stereo Pleiades have been acquired every year since 2013. These images are acquired in triplets and can be used to generate DEMs of resolutions better than 1 meter [Bagnardi et al. 2016]. However, these DEMs need to be carefully "cleaned" and post-processed to avoid artifacts that can greatly affect the modeling. Updating the DEMs is therefore done as often as possible, but it cannot be systematic.

Satellite sensor-derived TADR
The assessment of runout within our protocol depends on the reliability of the satellite-derived estimates of TADR, which are affected by uncertainties in the conversion between heat flux and volume flux as caused by thermal, rheological, and terrain factors [Harris and Baloga 2009]. Other uncertainties result from errors in the measurement of the thermal flux, mainly due to the presence of clouds that can attenuate the thermal signal detected by the satellite [Coppola et al. 2013], as well as pixel overlap and the point spread function at high scan angles [Harris et al. 1997;Harris 2013]. At Piton de la Fournaise, weather conditions are often cloudy, which therefore reduces the usability of satellite images to extract accurate TADR Thivet et al. 2020].
The method of discharge rate estimation by MIROVA differs from the approach of HOTVOLC. MIROVA system uses Low-Earth Orbiting platforms (Terra/Aqua) with pixels that are nominally 1 km 2 and requires a single image to provide a lava discharge rate. The error associated with cloud-free satellitederived time-averaged lava discharge rates with MIROVA is about ±50 % which can be reduced to ±35 % after calibration on previous eruptions ]. These TADRs represents an estimate of the volume of lava emitted in the hours preceding image acquisition, typically 6-18 hours in the case of data acquired in the mid-infrared [Coppola et al. 2010]. Consequently, the discharge rate obtained by MIROVA does not represent an instantaneous effusion rate value, but is a lava flux time-averaged over several hours [cf. Harris et al. 2007]. This temporal resolution (one image every 6 hours) does not allow for detailed evaluation of the lava discharge rate variation for short (a couple of days) eruptions as can happen at Piton de la Fournaise.
With HOTVOLC, discharge rate is estimated using lava volumes acquired every 15 minutes, as calculated from infrared data provided by the MSG-SEVIRI (Spinning Enhanced Visible and Infrared Imager) sensor [Labazuy et al. 2012;Gouhier et al. 2016;Gouhier 2022]. Using this geostationary sensor allows detailed time-series to be generated over short periods but the main limitation is the low spatial resolution, where pix-els are~24.5 km 2 at the image location of La Réunion and this prevents the detection of low amplitude thermal anomalies. HOTVOLC uses the pixel-integrated thermal anomaly measured at a given time that is the balance of contributions related to hot lava material newly emplaced and the cooling lava material previously emplaced [Gouhier et al. 2016;Dumont et al. 2021]. This method has uncertainties of the order of ±50 % and topographic shadowing [Dehn et al. 2002], as well as problems associated with radiance smearing due to the detector's point-spread function [Mannstein and Gesell 1991] and pixel overlap [Cahoon et al. 1992], are particularly troublesome at high scan angles like at Piton de la Fournaise [Harris 2013].
In our protocol, we use the two systems together for crosschecking and on-the-fly validation. This strategy is also adopted by LAV@HAZARD [Vicari et al. 2011;Ganci et al. 2012].
Future improvements include the development of ad-hoc satellite missions with a volcanological focus that would constitute a turning point in monitoring these types of events, allowing us to get closer to instantaneous values ]. In the meantime it is possible to increase the temporal frequency and spatial resolution of thermal images through multi-platform data analysis [Harris et al. 1998]. For example, the implementation of VIIRS (Visible Infrared Imaging Radiometer Suite) data within the MIROVA system will double the number of daily images useful for an accurate estimate of TADR [Campus et al. 2022]. Similarly, the analysis of high-spatial-resolution images, such as those from Sentinel 2 and Landsat 8 and 9 will increase the probability of obtaining maps useful for locating the vents and active lava fronts.

Timeline to map production
Running DOWNFLOWGO is a very rapid process. Once the operator of the model has the vent coordinates, the map production and delivery can thus be executed in about 10 minutes, as it was the case for the December 2020 eruption. In practice this time is often longer, because of readjustment of some parameters, the need for scaling, or addition of particular features if needed or requested by the director of OVPF or EMZP-COI. With the map template ready, the key factor to ensure an efficient hazard assessment is therefore the availability and accuracy on the vent location.
The time constraint for delivery of the map within the current protocol therefore mostly depends on the time between the start of the eruption and the confirmation of vent location. This may take from several minutes to a few hours, and depends on: 1) The location of the main vent or fissures, which may be difficult to confirm if located far from any viewpoints, roads or trails, or out of view of the monitoring camera network; 2) The weather conditions, where if it is too cloudy or foggy, the eruptive fissures may not be visible from the various viewpoints and a flight survey cannot be conducted; 3) The availability of the OVPF staff for carrying out field surveillance. The staff involves a small team Peltier et al. 2022], and staffing problems can and do arise if an eruption occurs in the middle of the night, at a weekend, during vacations or periods of illness or lockdown ; and 4) The helicopter used for airborne surveillance operations is operated by the Gendarmerie and also serves for mountain rescue services around the island. Thus, the aircraft may be unavailable in a timely fashion if already out on, or called for, a rescue mission which is prioritized over surveillance duties.
These time constraints may be significantly reduced if the eruptive site is, instead, located using the tremor maps that are produced every 15 minutes by the OVPF and which are available online . Inside the Enclos, the tremor map provides a relatively accurate vent location because this is where the monitoring network is the densest (there are 20 seismic stations within the Enclos within an area of 100 km 2 , Figure 1). The center of the tremor location is obtained by a triangulation method, which locates the source of the tremor, interpreted as the main fissure, with an error of hundreds of meters [Battaglia 2003]. Eruptions outside of the Enclos would be more difficult to precisely locate with the tremor map because of the scattered and less dense station distribution. Stations should be placed on the western flank of the volcano to improve triangulation in this area. Vent location may also be improved with the increase of satellite sensor spatial resolution that would enhance our ability to locate the vent and hence to quickly respond to an emergency crisis.

Crisis management for future eruptions at Piton de la Fournaise
Since the protocol has been in place and operational, the eruptions of Piton de la Fournaise have been characterized as "typical" with relatively low effusion rates that reach a maximum of 30 to 40 m 3 s −1 . All eruptions have occurred inside the Enclos where there is no population or infrastructure, aside from the belt road, hiking trials, and the OVPF monitoring network. For these reasons, the present protocol has been employed without need for adjustment for eruption location and source terms (see Section 6.2). However, as shown by the hypothetical scenario considered in Section 4.4, if a fissure opens outside of the Enclos to feed a high-effusion-rate flow (>100 m 3 s −1 ) into an inhabited area, the response protocol will become extremely important as a response-support tool and will need careful tracking to assess its performance. The time between the onset of the eruption and delivery of product will be also crucial in such an event, but we run into the circular problem as to how far we can trust runs based on "typical" source terms for an "atypical" scenario. However, we stress that the most important parameter is the location of the eruptive fissure and its timely provision because it quickly provides an initial approximation on the hazard. This first-approximation need may be pre-empted, to an extent, by precursor analyses such as locations of pre-eruptive seismic activity, deformation, and soil degassing anomalies (CO 2 ). Further improvements may also rely on simulations of, and expert elicitations for, hypothetical scenarios [e.g. Tadini et al. 2022].

Exporting the protocol and the approach to other volcanoes
This protocol could be exported to any other mafic effusive volcano with some modification, as follows. For export, as a first step the FLOWGO runout and TADR conversion routine needs to be calibrated and validated for the new site, although this has, to an extent, already be done in many cases and a well-tested methodology exists for model initialization [e.g. Harris and Rowland 2001;Harris et al. 2007;Harris 2015;Coppola et al. 2016;Ramsey et al. 2019]. As a second step, an up-to-date DEM is required to calibrate DOWNFLOW. This, also, already exists at a number effusive centers including Nyiragongo, Mont Cameroun (Cameroon), and Capo Verde [e.g. Favalli et al. 2009;Richter et al. 2016], and the methodology for model initialization is also well developed ]. As detailed here, once the protocol has been initialized, calibrated, validated, and refined so as to be fit-topurpose, then a rapid response service involving timely provision of lava inundation hazard maps, tailored to the volcano and eruption in hand, should be attainable. Ultimately and importantly, engagement with end-users (local authorities) is required to refine any resulting maps to the local area and agency needs; this is the best way to ensure that the shortterm hazard maps are well understood and are useful for preparedness and reducing potential risks.

CONCLUSIONS
Near real-time lava flow forecasting and monitoring is a necessary step in mitigating the impact to communities on the flanks of effusive volcanoes. We present here an operational protocol that has become an indispensable tool for crisis management at Piton de la Fournaise since its first full implementation in 2019. To illustrate this effective near real-time protocol, we present the response to the 25-27 October, 2019, 10-16 February, 2020, and 7-8 December, 2020, eruptions as well as to a hypothetical eruption outside the caldera. The successful development of this protocol was only possible by collaboration and open discussion between scientists and civil protection. The developed protocol and user consultation allowed production of a scientifically valid short-term hazard map that provides probable lava flow paths and runout distances, dependent on the lava discharge rate, in a format easily understood by the responding agent and tailored to civil protection needs. The numerical lava flow modeling used as part of this protocol can easily be revised, adapted, and improved using an update-to-date DEM and constrains of the flow thermorheological properties. Adapting this protocol to other effusive centers is feasible following a four-step process that involves: 1) Calibration of the topographic model for flow path and inundation area (in this case, DOWNFLOW); 2) Calibration of the thermo-rheological numerical model for flow runout (in this case, FLOWGO); 3) Calibration of satellite-data-derived discharge rate derivation (in this case, MIROVA and HOTVOLC), and 4) Establishment of a communication format and map design that suits local civil protection expectations and needs.
Such a lava flow hazard response protocol needs to produce deliverables that satisfy the needs and requirements of the agencies responsible for managing the effusive crisis, which can only be accessible through open discussion and involvement of all actors from the scientific base to the civil protection end-point.

AUTHOR CONTRIBUTIONS
MOC led this study and ran the lava flow simulation operational tool at each eruption. AH conceptualized the near realtime protocol. DC developed and communicated the MIROVA data. MG developed the HOTVOLC system. AP and SD contributed to the map designed. NV produced the lava flow outline. All authors contributed to implement the discussion and writing the article.