Lava flow hazard map of Piton de la Fournaise volcano

Abstract. Piton de la Fournaise, situated on La Réunion island
(France), is one of the most active hot spot basaltic shield volcanoes
worldwide, experiencing at least two eruptions per year since the
establishment of the volcanological observatory in 1979. Eruptions are typically
fissure-fed and form extensive lava flow fields. About 95 % of some
∼ 250 historical events (since the first confidently dated
eruption in 1708) have occurred inside an uninhabited horseshoe-shaped
caldera (hereafter referred to as the Enclos), which is open to the ocean on
its eastern side. Rarely (12 times since the 18th century), fissures
have opened outside of the Enclos, where housing units, population centers,
and infrastructure are at risk. In such a situation, lava flow hazard maps
are a useful way of visualizing lava flow inundation probabilities over
large areas. Here, we present the up-to-date lava flow hazard map for Piton
de la Fournaise based on (i) vent distribution, (ii) lava flow recurrence
times, (iii) statistics of lava flow lengths, and (iv) simulations of lava
flow paths using the DOWNFLOW stochastic numerical model. The map of the
entire volcano highlights the spatial distribution probability of future
lava flow invasion for the medium to long term (years to decades). It shows
that the most probable location for future lava flow is within the Enclos
(where there are areas with up to 12 % probability), a location visited
by more than 100 000 visitors every year. Outside of the Enclos,
probabilities reach 0.5 % along the active rift zones. Although lava flow
hazard occurrence in inhabited areas is deemed to be very low
(< 0.1 %), it may be underestimated as our study is only based on
post-18th century records and neglects older events. We also provide a
series of lava flow hazard maps inside the Enclos, computed on a
multi-temporal (i.e., regularly updated) topography. Although hazard
distribution remains broadly the same over time, some changes are noticed
throughout the analyzed periods due to improved digital elevation model (DEM) resolution, the high
frequency of eruptions that constantly modifies the topography, and
the lava flow dimensional characteristics and paths. The lava flow hazard
map for Piton de la Fournaise presented here is reliable and trustworthy for
long-term hazard assessment and land use planning and management. Specific
hazard maps for short-term hazard assessment (e.g., for responding to
volcanic crises) or considering the cycles of activity at the volcano and
different event scenarios (i.e., events fed by different combinations of
temporally evolving superficial and deep sources) are required for further
assessment of affected areas in the future – especially by atypical but
potentially extremely hazardous large-volume eruptions. At such an active
site, our method supports the need for regular updates of DEMs and
associated lava flow hazard maps if we are to be effective in keeping
up to date with mitigation of the associated risks.


as the caldera sensu stricto (Michon and Saint-Ange, 2008) and is named the Enclos Fouqué (EF) (Fig. 1b). It is a relatively flat plateau in the middle of which an asymmetric terminal shield has been built to reach an altitude of 2632 m a.s.l. and with slopes of between 15° (on the northern flank) and 30° (on the eastern flank). The top of the terminal shield is marked by a 1.1 100 ´ 0.8 km summit caldera ("Cratère Dolomieu") that last collapsed in 2007 (Michon et al., 2009;Staudacher et al., 2009;Derrien et al., 2020). To the west, Cratère Dolomieu coalesces with a smaller (0.38 ´ 0.23 km) crater ("Cratère Bory"). The eastern part of the Enclos, i.e., below 1800 m a.s.l., is divided into two areas (most probably formed by successive landslides ;Michon and Saint-Ange, 2008) where the slopes are very steep (>30°), hereafter named "Grandes Pentes" (French for "steep slopes"), and the flatter coastal area, hereafter named "Grand-Brûlé" (French for "widely-burned"; Fig. 1b). 105

Effusive activity and risk
As observed on most basaltic shield volcanoes (Dvorak et al., 1983;Tilling and Dvorak, 1993;Walker, 1988), eruptions at Piton de la Fournaise are produced by dyke/sill propagation following preferential paths leading to the concentration of eruptive fissures within, tangential-to or radially-around the summit crater and along three main rift zones (Bachèlery, 1981). 110 At Piton de la Fournaise, there are three rift zones, the southeast rift zone (SERZ), the northeast rift zone (NERZ) and the north 120° rift zone (N120), and these are the zones where effusive activity is mainly concentrated (Fig.1b). Two other zones with a high concentration of eruptive fissures have also been identified, these being the South Volcanic Zone (SVZ) and the Puy Raymond Volcanic Alignment (PRVA; Michon et al. 2015;Fig. 1b).
Eruptions are mainly effusive and expressed fed by en-echelon fissure sets, with each fissure being a few meters to a few 115 hundred meters long. Activity at the fissures tends to be Hawaiian to Strombolian in style, to feed lava flows that can extend all the way to the ocean, a distance of up to 12 km downslope (cf. Sect. 6). The average volume of lava emitted per eruption between 1970 and 2007 was about 10 ´ 10 6 m 3 . In 2007, activity was punctuated by a low-elevation flank (590 m a. s. l.) event in the Enclos, during which 140 to 240 ´ 10 6 m 3 of lava was erupted in 30 days (Roult et al., 2012;Staudacher et al., 2009). Such atypical "low flank" events can be of high risk (Lormand et al., 2020), and historically the 2007 120 event has been, volumetrically, the largest single effusive event witnessed at Piton de la Fournaise. Through April 2020, it has been followed by a further 28 more "typical" eruptions, each with a typical volume of 5 ´ 10 6 m 3 and being fed from vents at higher altitudes (OVPF reports ISSN 2610-5101). Lava flows can be channel fed (e.g., Harris et al., 2016;Rhéty et al., 2017;Soldati et al., 2018), to feed extensive compound lava flow fields within the Enclos (Fig. 1b), although tubes can form during longer-lived eruptions, as in 2007 . Lava composition is usually transitional, ranging from aphyric 125 basalt to oceanite and midalkaline basalt (Albarède et al., 1997;Lénat et al., 2012), with eruption temperatures in the range 1150-1190 °C (e.g., Boivin and Bachèlery, 2009;Harris et al., 2019;Di Muro et al., 2014). Lavas thus have a relatively low viscosity of the order of 10 2 -10 4 Pa s upon eruption (Harris et al., 2016;Kolzenburg et al., 2018;Rhéty et al., 2017).
Most effusive eruptions at Piton de la Fournaise occur within the Enclos. This is part of La Réunion National Park and is uninhabited. However, being a major tourist attraction, the Enclos does receive about one hundred thousand visitors per year. As a result the hiking trails in, and access roads to, the Enclos can receive heavy pedestrian and vehicular traffic (e.g., 129 000 hikers accessed the summit in 2011 according to Derrien et al., 2018). The Enclos also hosts the national road (RN2) which crosses the Enclos from north to south at a distance of 800 m inland from the coast and at an altitude of 130-70 m a.s.l. (Fig. 1b). This is the island belt road and, if cut, severely impedes communications and travel between communities in the south and north of the island (Harris and Villeneuve, 2018). Based on an analyses of the records available since the first 135 observation in 1640, Villeneuve and Bachèlery (2006) estimated that 95 % of the eruptions have occurred in the Enclos Fouqué caldera (i.e., above 1800 m a.s.l.). These are, hereafter, termed "proximal" eruptions as they are proximal to the terminal shield.
Fissures and associated lava flows produced by these proximal eruptions may cut hiking trails and generate a risk for the visitors (Derrien et al., 2018), and necessitates evacuation and closure of the Enclos (Peltier et al., 2020). Fissures may also open ion the Grandes Pentes or in the Grand-Brûlé, i.e., below 1800 m a.s.l.. These are termed "distal" eruptions and these can 140 cut the RN2, as has been the case for eight eruptions since 1979 (Fig. 1b). Eruptions outside of the Enclos, termed Hors Enclos eruptions -"hors" being French for outside, are high risk as they can open in or above inhabited areas near the coast, such as above the villages of Sainte-Rose and Saint-Philippe (Fig. 1b). They can also open in the inhabited highlands of the volcano, where there are the towns of La Plaine des Cafres and Plaines des Palmistes (Fig. 1b). For example, the eruption of April 1977, which was the first Hors Enclos eruption since 1800, caused evacuation of Piton Sainte-Rose, burnt down or damaged parts 145 of 26 buildings (including houses, a church, the police station and a gas station) and buried part of the Sainte-Rose municipality (Kieffer et al., 1977). A second Hors Enclos eruption in March 1986 caused around 400,000 euros of damage to houses, household contents, agriculture, roads and utilities in the municipality of Sainte Phillippe (Bertile, 1987;Morin, 2012).

Cycles of effusive activity 150
At Piton de la Fournaise, the location and frequency of effusive activity is controlled by the shallow (<2.5 km below the summit) stress field. According to Peltier et al. (2009), Got et al. (2013 and Derrien (2019), the evolution of the stress field at Piton de la Fournaise can follow cycles that have typical periods of one year to slightly more than a decade. These cycles typically start with one or a few eruptions at the summit, followed by one or a few proximal eruptions (with vents opening on the terminal shield or at its base). They may end with a large-volume, distal eruption at a distance of more than 4 km from the 155 summit, as was the case in 2007, or an Hors Enclos eruption, as was the case for the cycle that ended in 1977. In addition to these superficial cycles, lava composition seems to be affected by cyclic changes in the "fertility" of the source magma as observed since 1930, the date since which the petrology of the volcanic products has been well documented (Vlastélic et al., 2018). These deep source-related cycles are typically of longer duration (20-30 years) and are characterized by a continuous increase followed by a continuous decrease in source fertility, until a minima is reached which marks the end of a cycle and 160 the beginning of a new one (Vlastélic et al., 2018). High-volume distal eruptions associated to crater collapse occurring in 1931-35, 1961, 1986 and 2007 mean that three-quarters of such cycles are characterized by such "atypical", low flank events, where they typically occur just after the peak of source fertility (Derrien, 2019). This cyclic activity shows that eruptive fissures are more likely to open at lower altitudes and outside of the Enclos towards the end of superficial cycles, and particularly voluminous effusive events may be expected just after the peak in a deep-seated cycle. 165

Digital Elevation Models
Given the sensitivity of lava flow paths to topography (Favalli et al., 2005), for a target where the topography changes so often as is the case at Piton de la Fournaise, using a frequently updated, and indeed up-to-date, Digital Elevation Model (DEM) is essential. For this study we could use to three DEMs. The first was obtained in 1997 for the whole island as derived from SIR-C images. This DEM has a vertical and horizontal resolution of 0.15 m and 25 m, respectively, and was available on a local 170 projection (Gauss-Laborde Réunion) where elevations were based on the ellipsoid (Villeneuve, 2000). The second DEM was acquired in 2008-2009 by airborne Light Detection And Ranging (LIDAR) and has a 5 m resolution over the whole island, and 1 m resolution near the coast. The vertical resolution is 0.05 to 0.1 m (Arab-Sedze et al., 2014). Finally, the most recent DEM was obtained in April 2016 from optical Pléiades satellite images acquired in Stereo Triplet mode and processed by the French Centre National d'Etudes Spatiales using their S2P restitution code (de Franchis et al., 2014). This DEM is restricted 175 to the Enclos, which is where all topographic changes took place between 2009 and 2016, and was sampled at 5 m resolution for this study.

Lava flow recurrence time
As a first step in building our lava flow hazard map, we needed to prepare a complete and statistically robust eruption inventory, to estimate lava flow recurrence time on the volcano flank. The inventory of all eruptions since record allows 180 (beginning of the 18 th century since the first confidently dated eruption in 1708) has been detailed by several studies (Michon et al., 2013;Morandi et al., 2016;Villeneuve and Bachèlery, 2006). By compiling these studies and the exhaustive record since establishment of the OVPF Roult et al., 2012;Vlastélic et al., 2018), as well as recent OVPF reports (ISSN 2610-5101), we estimate that Piton de la Fournaise had about 250 eruptive events over a period of three centuries. About 95 % of these have occurred inside the Enclos (Table S1). 185 To produce the lava flow hazard map, we need to account for recurrence of lava flows and not the number of eruptive events. Moreover, an eruptive event may produce several lava flows that can be located several kilometers apart. Therefore, an individual flow is counted when it is either temporally or spatially separated from another flow. Therefore, an individual flow is counted and entered into the data base when it is either temporally or spatially separated from another flow. For any eruption at Piton de la Fournaise a fissure typically opens perpendicular to the slope, the magma erupts uniformly along the 190 fissure to feed several lava flow units simultaneously to form a flow field of many overlapping flow units. In such a case we count only a single flow, the longest, and do not consider all units that comprise the compound lava flow field in the data base (e.g., Walker, 1973). In the case where there is a pause in the eruption and new flows are emitted in a second phase of activity, https://doi.org/10.5194/nhess-2020-394 Preprint. Discussion started: 11 December 2020 c Author(s) 2020. CC BY 4.0 License. then the longest flow after the hiatus is also counted, so that the eruption has two flows associated with it in the data base.
Finally, if an eruption feeds lava flows at multiple, spatially distinct locations, each lava flow site is counted as a separate unit. 195 Following this counting strategy, since the creation of OVPF (late-1979) until the end of 2019, a period when the volcano has been continuously monitored so that the inventory can be trusted to be 100 % complete, the volcano erupted 77 times within which we can identify 128 distinct lava flows (Table S1). This translates to around 1.7 lava flows per eruption.
Given the fact that lava flows are more frequent in the Enclos than beyond its limit, we had to count the number of lava flows inside the Enclos based on a different period than for the flows outside of the Enclos. For this, we divided the 200 volcano into five regions: the Enclos, and the Hors Enclos, which was sub-divided into northeast flank (NE), southeast flank (SE) and the highlands along the N120 rift zone (N120; Fig. 2a) and the rest. Inside the Enclos, we consider a time period since 1931 (the date since which a rigorous record and mapping has been possible; Derrien, 2019). Until the end of 2019, this involved a total of 195 individual lava flows ( Fig. 2a; Table 1; Table S1). The recurrence time of lava flows within the Enclos is therefore estimated at one every ~5.4 months for the 1931-2019 period. Over this period, only three eruptions occurred 205 outside of the Enclos, and six lava flows were counted for these three events (three in 1977, two in 1986 and one in 1998). To calculate the recurrence time of Hors Enclos lava flows, six is a rather small number if we are to ensure a good statistical representation. We therefore increased the time period to extend back to 1708. Over the 1708-2019 period, fifteen lava flows were registered outside the Enclos ( Fig. 2b; Table 1; Tabler S1). Nine lava flows occurred on the southeast flank of the volcano.
Of these, five were observed (in 1774, 1776, 1800 and two in 1986). A sixth event was dated at 80 ± 35 BP using C 14 dating 210 on the carbon in the soil below the Piton Raymond lava flow (Vergniolle and Bachèlery, 1982), and three others were dated at 1726, 1765 and 1823 by measuring the base diameter of pioneer trees (Albert et al., 2020). Five other flows were observed on the northeast flank (one in 1708, three in 1977 and one in 1998). Within the less active N120 rift zone, the number of historical eruptions is not large. According to Morandi et al., (2016), eleven tephra or lava flow deposits can be dated since 2920±30 BP, but only one lava flow has been dated after 1708 ( Fig. 2b; Table 1; Table S1). This eruption (named Piton Rampe 14) was 215 dated by C 14 at 140±90 BP (Vergniolle and Bachèlery, 1982;Morandi et al., 2016). However, this eruption was not observed and has not been further studied. Another recent eruption is also suggested from some poorly characterized deposits (lapilli) close to the Trous Blancs area at 145±30 BP (Morandi et al., 2016), but no lava flow has been identified associated with this event. However, if we consider only one lava flow since 140±90 BP, or the eleven eruptions since 2920±30 BP, the eruption recurrence time does not vary significantly (one every 263 years and 209 years for the two periods, respectively). We therefore 220 assume just one eruption over our 311 year period of records for the N120 rift zone.
Overall, the inventory of Table S1 represents a minimum bound on eruptive activity because it is possible that other eruptions and lava flows were not observed nor have yet been identified in the geological record or dated. The minimum recurrence of Hors Enclos lava flows since 1708 is therefore estimated at one every ~21 years (Table 1). The relative occurrence probability is thus calculated by normalizing the number of lava flows per year within each region over the given 225 period. The resulting relative probability of occurrence for a lava flow is 97.99 % inside the Enclos and 2.01 % outside of the Enclos. Beyond the Enclos relative probability can be divided into 1.29, 0.67 and 0.20 % for the southeast, northeast and N120 rift zones, respectively (Table 1).
We here follow the same methodology developed since (Favalli et al., 2009c) and assume that future vents are more likely to open in areas where previous vents cluster. The probability of future vent opening is therefore determined on the basis of location of historical vents and scaled to the lava flow recurrence probability within each region (cf. sect. 3), under the general assumption that the characteristics of future eruptions will be similar to those of past eruptions. 235 Our inventory of vents for Piton de la Fournaise is based on the mapped scoria cone distribution (Davoine and Saint-Marc, 2016; Michon et al., 2015;Di Muro et al., 2012;Villeneuve and Bachèlery, 2006) and is here updated with all new vents formed between 2015 and the end of 2019 (Fig. 3a). Finally, any vents that were missing from this data base but were apparent from the mapping of Derrien (2019) were added. The vent density distribution (number of vents per unit area) was then obtained by applying a Gaussian smoothing kernel to the map of vent locations (Bowman and Azzalini, 2003;Favalli et al., 240 2012;Richter et al., 2016), with a bandwidth that is a function of the local vent density. Within the Enclos, the highest concentration is located to the southeast of the summit crater and is up to 51 vents per km 2 . To the east, the concentration decreases from 8 vents per km 2 in the Enclos Fouqué caldera, to 3 and 0.5 vents per km 2 in the Grandes Pentes, and the Grand-Brûlé, respectively. Outside of the Enclos, resulting vent densities range from 0.01 to 8 vents per km 2 in the rift zones, with a highest concentration of up to 21 vents per km 2 within the PRVA (as already noted by Michon et al., 2015). 245 We find that, of the 726 mapped vents, 45 % (327) are within the Enclos, while 55 % (399) are outside of the Enclos.
For the Hors Enclos vents, 19.6 % (142) are on the N120 rift zone, 13.5 % (98) are on the southeast flank, 6.7 % (49) are on the northeast flank, and 15.1 % (110) are elsewhere on the volcano flanks (Table 2). These distributions of vents per region differ significantly from the lava flow relative occurrence probability per region (Table 1). For example, while the N120 rift zone accounts for 19.6 % of the total number of vents (Table 2), only 0.2 % of the lava flows recorded since 1708 have occurred 250 in this region (Table 1). Instead, within the Enclos, where ~98 % of the lava flows have been emplaced (Table 1), we count only 45 % of the total number of vents (Table 2). This difference is mainly due to the period of time required for cones located in the various regions to disappear. In the Enclos, resurfacing processes and landscape changes, such as fissure opening, building of new scoria cones over old ones and burial by lava flows occurs at much higher rates than beyond the Enclos.
Therefore, the lifespan of scoria cones (vents) within the Enclos is considerably shorter than compared to that of scoria cones 255 outside of the caldera, where 65.8 % of the Enclos has been resurfaced at least once since 1931. Moreover, the selected time period (1708-2019) does not cover the full variability in eruption location and activity cycles documented by geological studies (Morandi et al., 2016 and references therein), where (in the geologic past) eruptions outside of the Enclos have, at times, been much more frequent than during the last 300 years. For this reason, the sum of probability density function for future vent opening within a given region of the volcano was set equal to the relative occurrence probability of lava flow in that region. 260 The resulting map of the probability density function of vent opening per unit area is given in Figure 3b. The highest probabilities exceed 10 %/km 2 with a maximum of 24 %/km 2 at the summit crater and across the proximal area of the terminal shield, while moderate values are obtained across the Grandes Pentes and in the Grand-Brûlé (0.01-0.5 %/km 2 and 0.003-0.01 %/km 2 , respectively, Fig. 3b). The northeast and southeast rift zones also have low to moderate values (0.5-0.003 %/km 2 ), while the N120 rift zone has a value of less than 0.003 %/km 2 . For the rest of the volcano the probability of vent opening is 265 very low, being less than 0.0001 %/km 2 . However, we note an area of relatively high values (up to 0.5 %/km 2 ) located at the PRVA (Fig. 3b).

DONWFLOW model calibration
Lava flows are gravitational flows that follow the steepest descent path defined by the underlying topography (Harris, 2013). 270 DOWNFLOW is a numerical code that computes the steepest descent paths while randomly applying a given vertical perturbation (Δh) at every pixel of the DEM (Favalli et al., 2005). By iteration over a thousands of runs (N), the code computes all the possible flow paths from the vent to the edge of the DEM through introducing random noise (with a value of ±Δh) into each run, so as to produce a probabilistic estimation of the lava flow inundation area. Thus, before applying DOWNFLOW to a given case, the key input parameters Δh and N need to be defined for the volcano and DEM in question (e.g., Favalli et al., 275 2012;Richter et al., 2016). Calibrating DOWNFLOW therefore consists of finding the parameters Δh and N that are able to best fit the model-generated and actual lava flow areas (Favalli et al., 2005). To do this at Piton de la Fournaise, a range of Δh (0 to 5 m) and N (100 to 10000) were applied to the first lava flows that occurred immediately after the date of DEM generation.
DOWNFLOW then computes the array of steepest descent paths out to the limit of the DEM, which, in our case, is the coast, but does not allow computation of lava flow lengths. Therefore, for the calibration, the simulations were cut at actual the length 280 of the flow under consideration (Fig. 4). Following Favalli et al., (2009b) and Richter et al. (2016), the best fit parameters can be found by comparing the actual, mapped, lava flow area (AR) with that generated by the DOWNFLOW simulation (AS): If this ratio (µ) is one then the two areas coincide perfectly and the simulation is valid. However, as the ratio trends towards zero the simulation becomes increasingly unrealistic. 285 We performed three calibrations, one for each of the three available DEMs (Figure 4). In total, we ran 70,000 simulations and found that input parameters are slightly different for each DEM.  (Walker, 1973). Alternatively, expected cooling-limited flow lengths can be calculated using a thermo-rheological model, such as FLOWGO (Harris and Rowland, 2001), that use theoretical flow cooling and crystallization properties to estimate the point at which forward motion is no longer rheologically possible (Rowland et al., 2004;Wright et al., 2005). Here, at Piton de la 300

Probabilistic lava flow hazard maps 310
To produce the hazard map for each DEM, a 100-m spaced grid of computational vents was defined for the whole of Piton de la Fournaise and a total of 126,000 simulations, of which 12,000 were inside the Enclos, were run using DOWNFLOW. From this database of simulations, it was possible to construct the lava flow hazard map which gives the probability, at each point, of inundation by lava flow upon the occurrence of the next effusive eruption (Rowland et al., 2005;Wright et al., 2008). The hazard at any pixel i is thus defined as the probability Hi that the pixel will be inundated by a lava flow, and is given by (Favalli 315 et al., 2012): where the sum extends over all possible vent locations j with the coordinates xj and yj, Δx and Δy define the pixel size and rVj is the probability density function of vent opening at location j (this is shown for the whole edifice in Fig. 3b). The mask of the simulations obtained with DOWNFLOW is given by Pij, where pixel i has the value one if inundated by a lava flow 320 originating from vent location j and has a value of zero otherwise (irrespective of the distance between j and i). Finally, PLij is the probability that a lava flow originating from vent j reaches pixel i along the calculated flow path, and is obtained by converting the frequency distribution of lava flow length for all flows up to relevant date (Fig. 5) into a probability.

Hazard map for the entire volcano
To be as up-to-date as possible, the hazard map for the entire volcano was derived using a combination of the 2016 DEM for 325 the Enclos area and the 2010 DEM for the rest of the volcano (Fig. 6). We can do this, because the 2010 DEM will not have been affected by topographic change outside of the Enclos. We run DONWFLOW using the best fit parameters as determined for the 2010-2019 period (Fig. 4c). In the Enclos, the lava flow length distribution was determined over 1931-2019, while for rest of the volcano, we considered the length of the flow units mapped by Bachèlery and Chevallier (1982) (cf. Sect. 7; Fig. 5,   Fig. S1). This is then convolved with the probability density function of vent opening from all the historical vents presented 330 in Figure 3b following Equation 2 (cf. Sect. 5).

Inside the Enclos
The resulting map (Fig.6) clearly shows that the highest probability of lava flow inundation at Piton de la Fournaise is located within the Enclos, where 47 % of the Enclos has a high (>1 %) probability of lava inundation. Within the summit crater, in 335 some places of the Enclos Fouqué and across parts of the Grandes Pentes, the probability is extremely high (>10 %). In the Grand-Brûlé area the probability is intermediate (up to 2 %) to low (0.1-0.5 %) due to the low probability of vent opening as well as the relative short lava flow length distribution (in the Enclos only 11 lava flows have reached the ocean since 1931).
We note that the there is a relatively high probability (of up to 2 %) that the flows will cut the belt road in some places.

Outside the Enclos
The highest probability of lava flow inundation outside of the Enclos is up to 0.5 % and is on the southeast and northeast flanks of the volcano as well as in the PRVA (Fig. 6). The less active N120 rift zone has a very low, but still non-negligible, probability of lava flow invasion (<0.01 %). Although the probability of vent opening is higher at greater altitudes (Fig. 3), the lava flow length is usually longer outside of the Enclos than it is inside (Fig. 5). This implies that the probability of lava flow inundation 345 remains the same as distance increases from the vent. This means that outside of the Enclos the distribution of the probability of lava flow inundation seems to depend mostly on the preferential flow path rather than the distance from the vent.

Evolution of hazard maps with DEM within the Enclos
Given the very high number of eruptions and the existence of three DEMs, we can compare the evolution of the lava flow hazard distribution within the Enclos with time and construction of new topography (Fig. 7). To do this we built three hazard 350 maps, one for each of the DEMs. That is, in 1997, 2010 and 2016. For each map, we take the appropriate lava flow length distribution, relative recurrence time and vent distribution from time of DEM generation back to 1931 ( Fig. 5; Table 3), with the DOWNFLOW calibration tailored to each case (Fig. 4). The percentage of vents in the summit craters versus those forming in the rest of the Enclos is roughly the same at the end of 2019 as in 1997 (that is around 3.4-6.1 % of the vents are within the crater, Table 3). This is important, because lava emitted from a vent opening in the summit crater will become entrapped in 355 the pit, and will not to be free to flow down slopes of the Enclos. The results also show that since 1931 the number of vents opening increased from 88 between 1931 and 2010, to 186 between 1931 and 2019 ( Table 3). The number of lava flows emitted also increased from 106 for the period 1931-1997, to 200 for the period 1931-2019. Instead, the relative probability of lava flow in the summit craters decreased from 27.9 % in 1997 to 24.1 % in 2019, reflecting that eruptions within the craters become rarer over the analyzed period, becoming non-existent since 2014. 360 The three hazard maps for the Enclos (Fig. 7) show that the highest probabilities (2 to >10 %) of lava flow inundation are concentrated in the Enclos Fouqué area (above 1800 m a.s.l.). Between the 1997 and 2010, we note the development of a high probability area to the north of the terminal shield, and a slight decrease in probabilities to the southeast. This is due to the high number of eruptions occurring to the north of the shield between 1998 and 2006 (see Fig. S2). On the 2016 hazard map, the southeast part of the Enlos Fouqué shows a smaller high probability area, but with higher values of up to 12 %. This 365 reflects a focusing of vent opening in this area over the period 2010-2016. Figure 7 also illustrates that the hazard maps were able to predict where subsequent lava flows actually occurred. Lava flows subsequent to each map were all emplaced in high probability areas, and most of the flows did not extend onto the moderate (5 %) probability zone. The exception is the 2007 eruption that occurred at a low altitude in the Grand-Brûlé, and which represents the largest eruption of the considered period. This is the only event to have occurred in a low (<1 %;) probability zone (Fig. 7). 370

Hazard map interpretation
According to Calder et al. (2015), uncertainties in hazard maps are mainly related to three issues. These are (Calder et al., 2015): "(i) the incompleteness and bias of the geological record and the extent to which it represents possible 375 future outcomes; (ii) the fact that analyses based on empirical models rely on a priori knowledge of the events; and (iii) the ability of complex computational models to adequately represent the full complexity of the natural phenomena".
Our hazard maps presented here in Figures 6 and 7 have been produced from available historical and geological records of when and where past lava inundation has occurred at Piton de la Fournaise, as well as vent location, lava flow location, 380 recurrence, and length. The quality and detail of these records improve with time, and are better inside the Enclos than outside.
Secondly, the maps are based on stochastic simulations of lava flow paths, the reliability of which depend on the quality, and up-to-dateness of our topographic model (DEM). This is an issue at a frequently active effusive center, where emplacement of new lava units causes the topography to be in a state of near-constant change. Here we thus discuss the uncertainties related to the interpretation of the presented maps and the extent to which they are adequate in assessing risk associated with future 385 eruptions.

Historical and geological records: representativeness of future eruptions
Our geographical distribution of vent density per unit area (Fig. 3b) is based on an inventory of 726 vents and includes prehistoric vents (scoria cones) identified in the geological record back to at least 57 ka BP (McDougall, 1971). The lava flow inventory and recurrence probability inside the Enclos covers a time span of 88 years (from 1931 to 2019), and outside the 390 Enclos the record spans 311 years (from 1708 to 2019). Given the frequency of activity, this allows a comprehensive data set of 220 lava flows to be used (Table S1). These relatively long periods and large numbers of cases, mean that the hazards maps produced here provide the future vent opening and lava flow inundation probabilities relevant at a temporal scale spanning decades to centuries, and are based on a statistically robust data set. They should therefore be reliable and trustworthy for long term land use planning and management. This includes agricultural practices, urban planning, road network planning, 395 assignment of protected areas, park-use planning, implementation of trail networks and installation of subterranean and above ground electric, sewage and gas networks; as well as targeting of potential zones where repair and replacement will be necessary (e.g. Tsang et al., 2020). Nonetheless, although, we argue that this hazard map is trustworthy in the long term, it may not be adapted for the short term, i.e., the next few years, over small spatial scales, i.e., hundreds of meters. Given the frequent resurfacing of the Enclos, the topography of the Enclos is highly variable in the short term, where emplacement of 400 new lava flow units can cause subsequent flows to advance down paths that are displaced 10s to 100s of meters over what would have been the case had the initial, pre-eruption, topography applied (Harris and Rowland, 2015). As shown in Harris et al. (2019), emplacement of new flow units, locally, can invalidate a DEM in a matter of hours, causing subsequent flow paths impossible to replicate. Thus, to accurately and precisely forecast exact lava flow paths the topography needs to be revised after each eruption, and even during long-lasting (weeks-to-months long) eruptions. As shown in Figure 6, the lava flow 405 inundation probability distribution changes over the time span of a few years. However, although precise hazard maps for such active centers need to be regularly updated, over the time span of decades to centuries, the topography remains broadly the same, where in the case considered here the Enclos has a slope profile that is downhill from east to west, guiding flow towards the coast (Fig. 6).
The interpretation of the hazard map may also need to vary with time as a system cycles through phases of more and 410 less frequent activity, and as the intensity of events within each cycle wax and wane. Activity at Piton de la Fournaise follows cycles that evolve over time scales of a few years to a decade or so (Derrien, 2019;Got et al., 2013;Peltier et al., 2009;Roult et al., 2012;Vlastélic et al., 2018). This cyclic eruptive pattern will control the location, intensity and magnitude of the eruptions, and hence the location, length and area of inundation of associated lava flows. Eruptive fissures tend to be closer to the summit craters (i.e., proximal eruptions are more probable) at the start of a cycle than at the end (when distal eruptions are 415 more likely to occur). Moreover, cycle durations will depend on the magma supply rate, where the lower the supply rate the longer the cycle (Derrien, 2019). Cycles also often terminate with a relatively long-lived, high volume, and high intensity effusive event (Coppola et al., 2017). Therefore, depending on the length of cycle and where, temporally, in the cycle we are, both the geographical distribution of the probability of vent opening and the length of associated flow will vary.
Piton de la Fournaise also experiences longer-duration (20-30 year long) cycles that are related to systematic evolution 420 in the fertility of the magma source (Vlastélic et al., 2018). Phases of increasing mantle fertility are related to increasing eruption frequency and erupted volume, as the mantle porosity increases thus promoting melt extraction. Such increasing phases of activity ultimately culminate in exceptionally voluminous, eruptions with high effusion rates, as well as summit collapse or pit-crater formation; as occurred in 1931in , 1961in , 1986in and 2007in (Vlastélic et al., 2018. Given that the lava flow length will be related to effusion rate and lava volume (Harris and Rowland, 2001;Malin, 1980;Pinkerton and Wilson, 1994;425 Walker, 1973), the lava length frequency distribution used to compute the hazard map is thus dependent on where, in terms of time, we are in such a longer-term source-related cycle. In addition, the collapse and pit crater formation events at the end of the cycle will strongly modify the topography in the proximal region.
Because our hazard maps are based on a database whose events are dominated by "typical" eruptions (i.e., only four or the 137 eruptions since 1931 are high volume source-related-cycle terminating events), they are representative only of such 430 typical cases and event scenarios. Therefore, the maps cannot be applied to assess the hazard associated with such "atypical" (source-related-cycle terminating) events. To illustrate this point, we take the 1997 hazard map (Fig. 7a). In this map, the 2007 lava flow is not well predicted where it occurred in a low (<0.5 %) probability zone. For the reasons mentioned above, the occurrence of atypical eruptions (such as the 1931, 1961, 1986 and 2007 events) cannot thus be foreseen by this particular set of hazard maps. This is a major issue in terms of risk because such atypical events have exceptionally high lava discharge rates 435 (which were sustained at more than 100 m 3 /s over 30 days in 2007; Staudacher et al., 2009), with flows advancing rapidly (within hours of the eruption onset) cut the belt road and enter the ocean (Harris and Villeneuve, 2018a). As a result, they are capable of rapidly inundating large areas to cause extensive damage. They thus represent the highest risk effusive event at Piton de la Fournaise, but are not accounted for in our hazard maps. Dedicated studies on the probability of occurrence of such high magnitude and intensity, but atypical, events needs to be conducted, and a separate set of hazard maps are required for 440 such events. Likewise, our analysis does not consider the poorly studied, but relatively recent (post-1708), long-lasting activity related to overflow from summit lava lakes, as was common between 1750 -1800 and around 1850 (Michon et al., 2013;Peltier et al., 2012). Our maps are, though, applicable to the most commonly effusive event scenario currently encountered at Piton de la Fournaise. However, they must be used and trusted with the above caveats in mind regarding the type of activity and effusive event to which they apply.

Numerical modelling
The DOWNFLOW stochastic approach computes the lava flow inundation area by summing N lava flow paths for a given random vertical perturbation (Δh) added to the DEM. As discussed above, the DEM thus needs to be regularly updated when a high frequency of lava flow emplacement events modifies the topography. The DOWNFLOW calibration requires a choice of two main parameters: the number of runs and the maximum range of Δh. The parameter space of Figure 4 shows that the 450 required number of runs (N) needed to achieve a good fit with the lava flow area (µ > 0.5) must, in this case, be at-least 2 000-5 000. Therefore, a safe choice for N to ensure statistically robust simulations is 10 000 in all cases, and this is used here to ensure model (and output map) robustness (Tarquini and Favalli, 2013). The calibration for Δh gives different results for preand post-2007 cases, being 5 m and 2 m, respectively. According to Favalli et al. (2005), 2Δh represents the characteristic vertical height of an obstacle that the flow can overcome. This implies that, since 2007, lava flow dimensions have changed: 455 they must have become thinner as they can now only overflow an obstacle if that obstacle is 3 m lower than prior to 2007.
Indeed, the lava flows in our data base are more voluminous, longer and thicker between 1997 and 2007 (mean length: 6480 m, mean volume: 49 ´ 10 6 m 3 , mean thickness 12.5 m) than between 2007 and 2019 (mean length: 2500 m, mean volume: 3.5 ´ 10 6 m 3 , mean thickness 4.3 m). This difference can be related to the cyclic activity. The period between 1997 and 2007 corresponds to end of a deep source cycle (Cycle 3 of Vlastélic et al., 2018). As a result, lavas were increasingly associated 460 with higher magma fertility and magma supply rates, to produce more voluminous flows at relatively high effusion rates.
Hence, following Walker (1973) and Malin (1980) units were longer and thicker, and thus able to overcome higher obstacles.
However, the period after 2007 relates to the beginning of a new deep source cycle (Cycle 4 of Vlastélic et al., 2018), with lower fertility and output rates. During this period the reverse is true: flows were less voluminous, and hence shorter and thinner and able to only overcome lower obstacles. The calibration parameters used for DOWNLOW, especially Δh, thus must 465 be determined for, and then applied to, a specific period of activity if the results are to be valid; and changed should there be an evolution in a cycle or eruption style.

Use for hazard mitigation at Piton de la Fournaise
In terms of hazard mitigation at Piton de la Fournaise, our maps are intended to provide information of value in planning actions implemented by local authorities and in building resilience to future eruptions, with a focus on protection of public 470 infrastructure. Figure 6 shows the locations of all roads and hiking trails that are the most likely to be covered by future lava flows, and the lengths of each that may need to be replaced.
For the island belt road, we see that 4.5 km of the 10 km that crosses the Enclos is in an intermediate probability (>0.5 %, Fig. 6) zone in terms of lava inundation. We also see that even if the probability is low, the probability of lava flow invasion is not negligible for several municipalities situated beyond the Enclos on the northeast and southeast flanks of the volcano (i.e., 475 Saint-Phillipe and Sainte-Rose, Fig. 6). In terms of hiking trails, although the main trail to the summit crater is in a low probability zone in terms of lava inundation, it is in a very high probability zone in terms of vent opening (>2 %/km 2 , Fig. 3).
Regularly since 2014, trails inside the Enclos Fouqué have been covered by lava, as for example in February 2019 when 400 m of trail very close to the viewing platform for the summit crater was buried (OVPF reports, ISSN 2610-5101).
In the case of an imminent eruption, our hazard map will thus complement the real time information that OVPF shares 480 with civil protection, who need to identify potential locations from where people may need to be evacuated and road sections to be closed (Peltier et al., 2020). In addition, land use plans (referred as Plan de Prévention des Risques in French law) at La Réunion do not currently take into account volcanic hazard. Our maps are thus also intended to aid and guide stakeholders in developing effective and comprehensive mitigation and land use plans, with the caveat that our maps are for a "typical" effusive event. Separate maps will need to be drafted in the future for atypical events, including high intensity events and prolonged 485 lava lake overflow at the summit, and this map should be updated whenever topography changes and a new DEM becomes available. However, we have here provided a flexible methodology that allows for ease of up-date.

Conclusion
Piton de la Fournaise is one of the most active effusive centers in the world, it being a volcano that has experienced more than one eruption a year since human settlement of the island at the end of the 17 th century, and two eruptions a year since the 490 volcano has been closely monitored by OVPF since 1979. The resulting database to calculate lava flow frequency and probability of future vent opening thus comprises 207 individual lava flows emplaced within the Enclos since 1931, 15 lava flows emplaced outside of the Enclos since 1708 and 726 vents over the entire volcano edifice. This large data base allows compilation of a statistically robust hazard assessment and production of probabilistic hazard maps for lava flow inundation.
Within the most active area of the volcano (i.e, the Enclos), a lava flow has been emitted every 5.3 months since 1931, 13 of 495 which have cut the island belt road, while beyond the caldera, an effusive event has occurred on average every 21 years over the last three centuries. Since 1931, two Hors Enclos eruptions since 1977 have caused property damage and twelve lava flows have cut the belt road.
Our lava flow hazard map uses these data, with lava flow path projections based on the stochastic model DOWNFLOW, to identify those areas that are most likely to be impacted by lava flows during any future eruption, and to 500 quantify the probability of inundation. Fundamentally we stress the need for frequent update of the DEMs at such frequently active volcanoes, where topography is constantly evolving to influence flow path, so as to produce a series of hazard maps that evolve with time and topography. In addition, given that the activity follows cycles that cause the location, erupted lava volume, effusion rate and flow dimensions to evolve with time, both the maps and the parameters used for modelling need to be applied on a case-dependent basis. We here produce maps that apply to effusive events that typify the first and second phases of a 505 cycle, but not the terminating high-volume events. Dedicated studies for probability occurrence of such high-volume events thus need to be conducted to complete the hazard assessment for all effusive scenarios. Our lava flow hazard map methodology is, though, intended as a flexible approach that can be applied to frequently active effusive centers, producing up-to-date maps that can be updated as eruption conditions evolve. This is essential support for informed land-use planning, as well as for use in crisis response and in drafting of hazard mitigation plans. 510

Data availability
All data shown in this study including the eruption inventory and GIS layers (georeferenced .shp files of the vents, fissures and lava flows as well as the calculated density function of vent opening and hazard maps) are freely available upon request at the OVPF or in the supplementary material. 520

Authors contribution
MOC lead this study, wrote the manuscript and made all figures and tables. MF executed the vent probability distribution, the DOWNFLOW calibration and computed the hazard maps; NV wrote the DEM section; AP, ADM and NV contributed to the eruption inventory; AD and NV draw the lava flow outlines; PB compiled the GIS data. All authors contributed to 525 implementing the discussion and to the writing of the article.