Validation of an integrated satellite-data-driven response to an effusive crisis: the April–May 2018 eruption of Piton de la Fournaise

Satellite‐based  surveillance  of  volcanic  hot  spots  and  plumes  can  be  coupled  with 17 modeling to allow ensemble‐based approaches to crisis response.  We complete benchmark tests on 18 an effusive crisis response protocol aimed at delivering product  for use  in tracking  lava  flows.   The 19 response  involves  integration of four models: MIROVA for discharge rate (TADR), the ASTER urgent 20 response  protocol  for  delivery  of  high‐spatial  resolution  satellite  data, DOWNFLOW  for  flow  path 21 projections, and PyFLOWGO  for  flow  run‐out.   We  test  the protocol using  the data  feed available 22 during  Piton  de  la  Fournaise’s  April–May  2018  eruption,  with  product  being  delivered  to  the 23 Observatoire du Piton de  la Fournaise via Google Drive.   The response was  initialized by an alert at 24 19:50Z on 27 April 2018.  Initially DOWNFLOW‐FLOWGO were run using TADRs typical of Piton de la 25 Fournaise, and revealed that flow at >120 m/s could reach the island belt road. The first TADR (10– 26 20 m/s) was available at 09:55Z on 28 April, and gave flow run‐outs of 1180–2510 m.   The  latency 27 between satellite overpass and TADR provision was 105 minutes, with the model result being posted 28 15 minutes later.  An InSAR image pair was completed six hours after the eruption began, and gave a 29 flow  length of 1.8  km;  validating  the  run‐out projection.   Thereafter,  run‐outs were updated with 30 each new TADR, and checked against flow  lengths reported from InSAR and ASTER mapping.   In all, 31 35 TADRs and 15 InSAR  image pairs were processed during the 35‐day‐long eruption, and 11 ASTER 32 images were delivered. 33


Introduction 34
Throughout the 1990's and 2000's methods were developed to extract lava flow discharge rates from 35 1 km spatial resolution satellite data collected by satellite sensors operating in the thermal infrared 36 (e.g., Harris et al., 1997;2007;Harris & Bologa 2009;Coppola et al., 2010). At the same time, high 37 spatial resolution (30 m) satellite data were shown to be of value for mapping lava flow fields (e.g., 38 Flynn et al., 1994;Wright et al., 2000;Lombardo et al., 2009), with InSAR data allowing estimation of 39 lava flow areas, thicknesses and, hence, volumes (e.g., Zebker et al., 1996;Rowland et al., 1999;Lu et 40 al., 2003). In parallel, a series of lava flow models were developed to allow flow inundation areas to 41 be simulated (e.g., Young & Wadge, 1990;Crisci et al., 2003;Vicari et al., 2007). Increasingly, the 42 capabilities have been merged to allow an ensemble-based approach whereby satellite data from 43 multiple wavelengths and spatial resolutions are combined to allow maximum constraint and cross-44 validation (e.g., Patrick et al., 2003;Rowland et al., 2003;Wright et al., 2005) and source term input 45 into real-time lava flow emplacement models (e.g., Wright et al., 2008;Vicari et al., 2011;Ganci et 46 al., 2016). Since 2015, just such a response model has been developed at Piton de la Fournaise 47 (Harris et al. 2017), where we here review and validate an updated version of the protocol so as to 48 review an ensemble approach to responding to an effusive crisis. 49 The response protocol is based on in situ observations and data acquisitions carried out routinely by 50 the Observatoire du Piton de la Fournaise (OVPF) team and the integration of four models: MIROVA 51 (Coppola et al. 2016), the ASTER (Advanced Spaceborne Thermal Emission Radiometer) urgent 52 response protocol (Ramsey, 2016), DOWNFLOW (Favalli et al., 2005) and FLOWGO (Harris and 53 Rowland, 2001). MIROVA is a near-real time hot spot detection system that uses MODIS data, and 54 has been calibrated for TADR calculation at Piton de la Fournaise by Coppola et al. (2010), the ASTER 55 urgent response protocol is a means of automatically prioritizing and targeting ASTER data 56 acquisition during a volcanic eruption. Instead, while DOWNFLOW is a stochastic model that 57 assesses potential flow paths based on iterative runs over a DEM with random noise added, FLOWGO 58 can calculate the cooling-limit of flow down each path (Rowland et al., 2005;Wright et al., 2008). To 59 estimate the maximum distance a flow can extend at a given effusion rate, FLOWGO tracks the 60 thermal and rheological evolution of a control volume of lava as it moves down a channel, tracking 61 the volume until the volume cools and crystallizes to such an extent that forward motion becomes 62 rheologically impossible (Harris and Rowland 2015). FLOWGO has been initialized for and tested for 63 lava channels at Piton de la Fournaise by Harris et al. (2016) and Rhéty et al. (2017), and-to allow 64 improved model initialization, iteration and application-has been rewritten and rebuilt in Python as 65 PyFLOWGO (Chevrel et al., 2018). It is this version of FLOWGO that we use here. 66 As described in Harris et al. (2017), the response protocol is initialized with the alert of an imminent 67 eruption and provision of the vent location provided by the OVPF as part of their mandated 68 monitoring and response procedures. Subsequently, it involves calling each model in sequence and 69 passing results between each actor, and then final product to OVPF, in as timely fashion as possible. 70 The protocol also calls in ground truth (for vent locations, effusion rates, channel dimensions, flow 71 lengths) provided by the OVPF as well as textural and chemical data (for eruption temperatures, 72 vesicularity, crystallinity, rheological models) produced at LMV, to improve model uncertainty and 73 syn-response validation. We show here how the response protocol works, and define the main 74 uncertainties, using a real-time exercise held immediately after the April-May 2018 eruption of Piton 75 de la Fournaise. The aim of the exercise was to refine model initialization and execution for Piton de 76 la Fournaise, reduce uncertainty, and to fully define the call-down and communication protocol. It 77 involved first following the data feed and executing responses, in the order that they were received, 78 followed by a validation phase in which remote sensing and model based estimates for discharge 79 rate and flow length were compared against ground truth. In doing so, we show how an integrated 80 multi-sensor remote sensing approach can be used to follow, document and quantify an effusive 81 event in near-real time. 82 The April-May 2018 eruption of Piton de la Fournaise and available data 83 The April-May 2018 eruption of Piton de la Fournaise began late on 27 April (19h50 UTC) from five 84 north-south orientated en-echelon fissures that opened between the elevations of 2165 m and 2285 85 m on the southwest flank of the terminal cone ( Figure 1a). Initially flow was channel-fed 'a'a which 86 moved down the SW flank of the Dolomieu. In a short time activity reached a peak and became 87 focused at a main vent roughly central to the fissure line at an elevation of 2200 m. Another much 88 less active vent a few meters to the north continued to project tephra and emit flames. Around the 89 two vents, scoria cones and tephra fields were constructed. Upon reaching the base of the Enclos 90 Fouqué wall (between the 29 and 30 April), lava flows turned southeast to follow the base of the wall 91 reaching a distance of 2.6 km before discharge rates declined and active flow fronts retreated to 92 positions closer to the vent (Figure 1b). Between 4 and 7 May, flow activity was concentrated in the 93 proximal section of the flow field with several tubes and, with two main zones of breakout being 94 active 200 and 500 m down the tube system ( Figure 1c). Breakouts from the tube system fed low-95 discharge rate flows which extended no more that 100-200 m. From 7 May new lava flows broke out 96 from an ephemeral vent at the base of the Enclos Fouqué Southern wall producing local vegetation 97 fires. Over the following days, the tube continued to extend and feed lava flows from its terminus, so 98 that by 10 May the tube exit was around 3.2 km from the vent. This continued to feed low-discharge 99 rate flows that extended over 1.1 km (or 4.5 km from the main vent) along the base of the Enclos 100 Fouqué wall. Activity continued in this way until 1 June 2018 when activity died out around 14h30 101 (local time). During the 34.6-day-long eruption, six aerial photograph, two aerial IR image and 102 several field observation campaigns, including GPS measurements, lava and tephra sampling, gas 103 analysis and UAV over flights were completed by the OVPF. In addition, 35 cloud-free MODIS images, 104 11 ASTER images and 15 InSAR image pairs all of which were available for near-real time analysis and 105 reporting. 106

Methodology 107
While implementation of MIROVA and the ASTER urgent response protocol (URP) allow near-real 108 time collection and processing of satellite thermal data for derivation of time-averaged discharge 109 rate and mapping of a thermal anomaly, DOWNFLOW and FLOWGO (DOWNFLOWGO) allow the flow 110 paths and potential run out distance to be projected. These models are called in sequence, where 111 the call-down procedure is given in Figure 2. As part of this system, output and product are shared 112 using a standardized reporting form (as given in Appendix A) which is shared between an email 113 distribution list involving all actors in the response chain, and to OVPF for integration into 114 surveillance and reporting duties. With each update, the group is issued an update email, flagging 115 the field that has been updated and giving the time and date of the update as well as the name of 116 the person responsible for the update. The reporting form has four fields for: (i) current MIROVA-117 derived TADR and time series; (ii) current vent location and DOWNFLOWGO projections; (iii) current 118 ASTER thermal distribution map, with flow field evolution time series and report; (iv) InSAR-based 119 flow length report and coherence images (Appendix A). Another field may be added to the reporting 120 form including relevant OVPF data collection, e.g., flow length from Structure from Motion (SfM), 121 SO2 flux, sampling locations etc. This is left at the observatory's discretion to add depending on work 122 loads and time commitment. 123 MIROVA and ASTER were called using the observatory bulletin announcing implementation of alert 124 level 1, that is an eruption is believed (on the basis of seismic and ground deformation data) to be 125 "imminent" (in the next minutes/hours). This causes ASTER to be targeted, and MIROVA to set up a 126 "watch" for the first sign of a hot spot. Upon eruption onset, DOWNFLOWGO is run as soon as vent 127 location(s) (GPS coordinates) is (are) known. The first vent location is usually provided by OVPF 128 personnel or gendarmerie using hand-held GPS from a helicopter which is flown by the police 129 (gendarmerie) service. Precision may vary depending on flight time available, the height of the 130 fountains and the number of aircraft in the air space above the eruption site. Initially, to give an 131 immediate idea of likely flow paths and inundation areas, 10000 flow lines are run to the edge of the 132 DEM (i.e., the coast) over the most recent 5-m DEM with random noise of between ±0.8 m and ±2.5 133 m being added between each run. The slope from the line of steepest descent (LoSD) at ±0.001 m is 134 then extracted (and smoothed every 10 m) and used for preliminary FLOWGO runs at various 135 effusion rates (10,20,30,40,50, to 100 m 3 /s). To do this, FLOWGO is initialized prior to the call 136 down using typical Piton de la Fournaise thermo-rheological conditions and textural properties as 137 given in Table 1. At the beginning of the eruption, a typical channel width of 4 m is taken ( Upon receipt of the first ASTER imagery a thermal anomaly map is produced, and flow locations and 146 lengths assessed on the basis of the spatial distribution of spectral radiance in 90 m ASTER band 12 147 (thermal infrared, 8.925-9.275 µm). In addition, vent location is checked where the intense thermal 148 anomaly at the vent is apparent in ASTER band 3 (near-infrared, 0.807 µm) image. The 15 m-pixel 149 size, and one pixel accuracy of the geolocation, allows the location of the vent hot spot to ±15 m. 150 This is often better than that provided by hand-held GPS, which when run in a fast moving helicopter 151 records a point that will lag behind the craft point by several hundred meters. If this is the case, the 152 vent location is updated and new DOWNFLOWGO runs are produced. If tubes begin to extend from 153 the vent, this-following Wright et al. (2000)-becomes apparent in the high spatial resolution 154 satellite images from the distribution of spectral radiance. In such as case, the source for 155 DOWNFLOWGO will be moved to the tube exit. 156 In addition, InSAR interferograms and SfM data are processed for flow thickness and length maps 157 that both add to the information flow and allow validation of model-based flow-length projections. 158 Although remaining largely underutilized in an operational response sense, the value of such data in 159 producing lava flow thickness maps as long been known (e.g., Zebker et al., 1996;Rowland et al., 160 1999;MacKay et al., 1998;Stevens, 2002;Lu et al., 2003), as has the potential for merging with 161 ancillary data, such as thermal-IR-derived TADRs and model-based lava flow run-outs (Rowland et al., 162 2003). The InSAR method consists of computing an interferogram by subtracting the phase between 163 two SAR images acquired for the same area at different times (for details of the method see 164 Appendix B). These statistics which are input into a fourth field in the reporting form (Appendix A) 165 and are also used to update the DEM used for flow path runs. 166

Validation 167
On 4 May 2018 an over flight was made in an ultra-light aircraft at a flight height of around 310 m 168 above the ground surface. A thermal camera was used to collect 52 images of the lava flow field and 169 vent system between 12:15 and 12:30 local time. The thermal camera was a FLIR Systems T650 170 which provides a 640 × 480 pixel image in the 8-14 µm waveband, with 0.65 mrad pixels. This, over 171 a line-of-sight distance of 460 m (and viewing angle of 48°) gives a pixel size of 0.3 m. Images were 172 used to obtain vent (eruption) temperatures and down channel surface temperature profiles to use 173 in FLOWGO, as well as channel and flow dimensions plus radiative (Q rad ), convective (Q conv ) and total 174 (Q tot = Q rad + Q conv ) heat fluxes to check against model output. In addition, the MODIS and ASTER 175 images collected at 10:30 (local time) on the same day (i.e., two hours previously) were fitted to the 176 thermal camera image mosaic to allow the heat fluxes and TADRs to be compared. TADR was 177 extracted from the thermal camera images using TADR = Q tot / (c p T + f ), in which is the lava 178 density, c p is specific heat capacity, T is the cooling range, f is the fraction of crystals grown down 179 flow and is latent heat of crystallization. Values characteristic of recent lavas at Piton de la 180 Fournaise were used for , c p , and f, these being 2079 kg/m 3 , 1225 J/kg K and 0.1, respectively, with a 181 cooling range of 75-250 °C (Harris et al. 2007). At the end of the eruption, following sample analysis, 182 the chemical, temperature, crystallinity and vesicularity sections of the initialization file for flow 183 modeling are checked, and if necessary, updated (Table 1). 184

Results 185
The trigger for the protocol of Figure 2  produced and posted on the reporting form ( Figure 3). This revealed that flows fed at sustained 207 rates in excess of 120 m 3 /s were capable of reaching the island belt road, to reach the coast. 208 However, because a 4 km wide basin existed after a distance of 4 km from the vent, flows became 209 held up at this point, with even flows at 80 m 3 /s coming to a halt 4 km from the vent; and to push the 210 model across the basin needed more than 120 m 3 /s. Thus, in reality, our prediction was that either 211 time would be needed to fill this basin, where lava needed time spread and pile up, and/or for a tube 212 to develop across the basin-a little like the case of lava flow advance towards Etnea Zafferana in 213 1992 (Barberi et al. 1993). 214 The first cloud-free MODIS overpass occurred at 09h55 (UTC, 13h55 local time) on 28 April, i.e., 215 around 14 hours after the eruption began. This yielded a TADR of 10-20 m 3 /s (Table 2). These 216 values were immediately input into the reporting sheet, thereby being handed onwards for input 217 into the PyFLOWGO initialization file. The first lava flow projection map was thus also completed and 218 posted; revealing flows were capable of extending up to 1180-2510 m under initial conditions 219 ( Figure 4a). The latency between satellite overpass and TADR provision was 105 minutes, with the 220 model result being posted 15 minutes later. The first S1B InSAR image pair was completed around six 221 hours after the eruption began and was also entered into the reporting sheet ( Figure 5a). These 222 revealed that the flow was already 1.8 km long and covered an area of 0.5±0.1 × 10 6 m 2 (Table 3);  223 giving an initial extension rate of around 5 m/min and coverage rate of 1400 ²/min. On the same day, 224 at 09h00 (local time), the first SfM survey was completed and by 16h00 (local time) approximate 225 location of the fissures and flow outline from aerial images were published by the OVPF. 226 At 03h33 (UTC, 07h33 local time) on 30 April, after a new aerial visit of the eruption, the center of 227 the main fissure was precisely given at 365365 m E; 7648810 m S. By this time, however, MODIS-228 derived TADR had declined to 3.7-6.9 m 3 /s (Table 2). Updating PyFLOWGO revealed reduced run-229 outs of 0.7-1.0 km. The first cloud-free ASTER image was acquired on 2 May. This revealed an 11 230 pixel-long anomaly of saturated pixels orientated NE-SW on the south flank of the Dolomieu ( Figure  231 6)-equivalent to a 990 m long zone of active lava (Table 4). The active vent was apparent as a single 232 pixel hot spot in the 15-m near-infrared data and the vent location was updated to 365280 m E; 233 7648835 m S (Appendix D), with the TADR for this day being 3.6-4.6 m 3 /s (Table 2). These details 234 were updated in the reporting form, and the vent location for DOWNFLOWGO adjusted slightly 235 (although this had no effect on the flow paths or LoSD). The following day (3 May), the second 236 coherence map was produced. This revealed that the lava flow field had, at some point, reached the 237 base of the caldera wall, turning SE to follow the base of the wall ( Figure 6) having attained a length 238 of 2.5 km ( Table 3). The shorter length of the active flows implied by the size of the thermal anomaly 239 in ASTER on 2 May, as well as the 17 pixel (1530 m) long zone of cooler pixels beyond the front of the 240 main hot spot indicated that flow front locations had begun to retreat back up flow by this time. 241 The thermal camera imagery obtained from the over-flight of 4 May confirmed that activity had 242 diminished, and comprised tube-fed breakouts of channel-fed 'a'a ( Figure 1c). Two main breakouts 243 were located where the southern breakout was fed by a 2 m-wide channel which fed a 110 m pad of 244 'a'a. Total heat flux from the breakout was 435±50 MW, which converted to a TADR of 0.61-1.65 245 m 3 /s. ASTER imagery revealed that, by 9 May, the tube had extended 2430 m ( Figure 6) to feed lava 246 flows of around 1.4 km in length. At the same time, MIROVA revealed continued decline in TADR 247 ( Figure 7) to between 1 and 2 m 3 /s. As a result, the vent location for DOWNFLOW was moved to the 248 tube exit, which ASTER gave as being at 364685 m E; 7647090 m S, and FLOWGO run at 1.6 and 3.8 249 m 3 /s (Table 2) with an updated channel width and eruption temperature (Table 1). This gave flow 250 lengths of 1-2 km beyond the end of the tube system ( Figure 4b). 251 Thereafter, TADRs remained at low levels ( Figure 7) and the flow field continued to build parallel with 252 the base of the caldera wall ( Figure 5). TADRs of 0.8 m 3 /s characteristic of the final week of the 253 eruption (Table 2) gave flow lengths that extended just 1 km from the end of the tube (Figure 4b). 254 The flow field (both predicted and measured) attained a final length of 4.1 km and area of 1.3±0.1 × 255 10 6 m 2 (Table 3), and a volume of 5.5±1.6 x 10 6 m 2 (Table 2). In all, 35 TADR sets were processed by 256 MIROVA (Table 2), 15 image pairs were processed for coherence analysis (Table 3), 11 ASTER images 257 were obtained using the ASTER URP (Table 4) and DOWNFLOWGO was launched three times as TADR 258 and vent location changed. Additionally, six aerial photograph data sets, two aerial IR image surveys 259 and multiple field observations, including lava and tephra sampling, gas analyses, UAV over flights 260 were completed by the OVPF. The final reporting sheet, filled out with all derived values from this 261 data set, is given in Appendix D. 262

Discussion 263
The aerial survey mapping of the flow field of 30 April allowed checking of the dimensions of the lava 264 flow field derived from InSAR data; the center line length being 1.8 km (the same as that given by 265 InSAR) and the area having excellent coincidence with the zone of incoherence obtained from the 266 InSAR data. Likewise, dimensions of InSAR zones of incoherence, ASTER thermal anomalies and 267 FLOWGO lengths are in excellent agreement ( Figure 8). For example, the thermal anomaly in ASTER 268 on 2 May revealed that flows had extended to a maximum distance of 2520 m in the proceeding 269 days. This compares with the 2.5 km long zone of incoherence recorded by the InSAR pair processed 270 the following day (3 May) and the 2510 m flow length generated by FLOWGO using the maximum 271 TADRs obtained from MIROVA the first few days of the eruption. Closing the circle with validation of 272 the FLOWGO run outs with good fits with dimensions of incoherence and thermal anomalies in InSAR 273 and ASTER data gives us confidence in the source terms (including MIROVA-derived TADR) entered 274 into the model. We next assess the uncertainty in those MIROVA-derived TADR, as well as the 275 FLOWGO run-outs and errors due to DEM problems. 276

Validation of MIROVA-derived TADR 277
The image collected on 4 May by MODIS-MIROVA indicates a total radiant power (Q rad ) of 497±149 278 MW, corresponding to a total TADR of 2.6±0.6 m 3 /s. Total radiant power is around 42 % of that 279 measured for the south breakout on 4 May using the thermal camera (i.e., 209±20 MW). The TADR 280 (1.13±0.52 m 3 /s) obtained from the thermal image is also 43 % that of the MODIS-MIROVA, 281 indicating confidence in the latter value and the conversion routine used. The third uncertainty is on surface temperature which controls heat loss and hence cooling rate. We 309 have used the typical effective radiation temperature approximated from the data of Flynn and 310 Mouginis-Mark (1994) for a lava channel on Kilauea to initialize the model ( Our final uncertainty is on channel width. If, for example, we reduce to a width of 2 m, depth and 316 velocity have to increase to 1.1 m and 4.8 m/s to balance the TADR. This yields a runout of 2550 m 317 or 46 % longer, so that an uncertainty on channel depth of 50 % yields uncertainty on runout of the 318 same order. However, to extent uncertainties may cancel. If for example, we increase the 319 vesicularity to 50 vol.% but decrease the phenocryst content to 1 vol.% we change the runout by just 320 50 m (for the same TADR). Likewise, if we decrease channel width to 2 m, but increase surface 321 temperature to 740 °C we change the runout by 50 m. Thus, our error appears to be around 4-5 %, 322 so that the error on a predicted runout of 3000 m, is just less than a few hundred meters. 323

DEM uncertainty 324
Until now, for the near real time response at the effusive crises at Piton de la Fournaise, DOWNFLOW 325 was run on the SRTM DEM from 2005. When we first ran the DOWNFLOW simulation (in May), the 326 LoSD was not south to the base of the Enclos Fouqué wall, but projected due East. It was not possible 327 to simulate the actual flow path because post-2005 topography could not be accounted for. 328 However, now that we have updated our flow projection by using the 5-m 2010 DEM to which lava 329 flow fields from October 2010 and August 2015 (which were both in the southern area of the Enclos 330 Fouqué) were added, the predicted path is south, moving around the western edge of the 2015 flow 331 field, and to reach the wall before flowing to the east along its base. This was exactly the trajectory of 332 the flow. Note that although the eruptive fissures were located near and onto the February 2015 333 flow (on the distal part) we did not update the DEM with this lava flow as it did not interfere with the 334 ongoing flow process. To model flows on a very active volcano, where topography is constantly 335 changing, we thus need a DEM that is updated after each eruption, so as to reduce uncertainties on 336 predicted flow inundation area. 337 To obtain the inundation area, DOWNFLOW needs to be calibrated to a specific scenario, and this is 338 achieved by tuning N and Dh (Favalli et al. 2011 In the present case, the changes in vent location between the first estimation and the coordinates 352 obtained in the field or from the satellite images did not change significantly. The effect on the 353 predicted flow path was therefore minor and limited to within 100 m of the vent area. However, 354 knowing, and moving to, the break out location of 9 May, was essential to predict the final flow 355 length (at the given new TADR). The protocol we are offering here, that is sharing ASTER, MODIS, and 356 DOWNFLOWGO allow a back and forth to update the vent location and is therefore of major 357 improvement for correct estimation of the lava flow path and runout distance. In addition this 358 protocol is of service to OVPF to aid in monitoring needs for lava flow field evolution allowing both 359 crisis management and appraisal of need to evacuate ground based monitoring stations falling in 360 flow paths. 361

Conclusion 362
With the near-real time availability of data from so many satellite-based sensors, as well as the 363 immediate availability of ground truth through upload to internet-based data hubs, the best way 364 forward to tracking an effusive crisis is an ensemble-based approach. we have begun to convolve data from sensors flown on UAVs, as well as from crowdsourcing. 371 Another developing avenue is small, low cost satellite networks, such as the small satellite 372 Technology Experiment Carrier-1 (TET-1) as developed by the German Aerospace Center and 373 dedicated to monitoring high temperature events (Zhukov et al. 2006). Such systems offer high 374 spatial resolution (160 m) thermal infrared imagery at a relatively high temporal resolution (3 days) 375 and have shown to be of value in tracking effusive crises yielding TADR time series to supplement 376 those provided by MODIS (Zakšek et al. 2005). 377 As shown here, merging thermal data of different resolutions allows time series generation with the 378 best possible temporal resolution and precision; cross-validation of TADR, error and uncertainty 379 assessment; and input into lava flow emplacement models. The next step will be the use of InSAR 380 data to allow DEMs to be updated between eruptions so as to ensure that flow paths are correct and 381 use the most up-to-date topography available, with the DEM evolving as the topography changes. 382 This is a key feature, especially during a long-term eruption with changing topography and vent 383 position. In turn, the chain can be inverted where good agreement of model-predicted flow lengths 384 with dimensions of thermal and incoherence anomalies in high spatial resolution and thermal data 385 suggests that the source terms input into the model are valid.  Yellow stars give the distance down the LoSD FLOWGO runs at each generic effusion rate 523 (numbers are in m 3 /s). These are the "effusion rate contours" for this eruption. 524 InSAR incoherence (blue outline) and field mapping (yellow outline) on the same dates. 528 Background shows the DOWNFLOW inundation area. 529    Table 1. Key thermal, textural and rheological source terms used to initialize PyFLOWGO at Piton de la Fournaise as given by Chevrel et al. (2018). These are based on measurements and best-fit testing of FLOWGO on lava channels active during the December 2010 eruption of Piton de la Fournaise as described in Harris et al. (2016).

InSAR data processing
The InSAR method consists of computing an interferogram by subtracting the phase between two SAR images acquired of the same area at different epochs. The phase recorded on a SAR image depends both on the radar wave round trip travel time between the instrument and the ground, and on the interaction between the radar wave and reflectors on the ground surface. Provided that this last component remains stable between two successive acquisitions, the differential phase displayed on the interferogram will only reflect changes in the radar wave travel time between the two acquisitions and can be exploited to measure possible ground surface displacement occurring between the two acquisitions. This is the classical application of InSAR that has seen a succesful development in the field of volcanological application since the pioneering work of Massonnet et al. (1995) (e.g., Biggs andPritchard 2017;Pinel, Poland, and Hooper 2014;Massonnet, Briole, and Arnaud 1995;Ferretti, Prati, and Rocca 2001;Hooper et al. 2004).
If the geometry or dielectrical properties of the reflectors on the ground change significantly between the two radar acquisitions, the differential phase will appear very noisy ("decorrelated" or "incoherent") on the interferogram, making it unexploitable for displacement measurement. This generally occurs when the ground is covered by dense vegetation, the geometric properties of individual plants changing very quickly due to continuous growth or to the motion of leavees in the wind. It could also occur when heavy rain, snow fall, strong erosional events, air fall (in volcanic context) occurs between the two radar acquisitions. Interferometric coherence provides an estimation of the temporal stability of the ground contribution to phase measurement. It is generaly derived as an inverse function of the phase variance calculated between neighboring pixels (e.g., in a 3 × 3 pixels box) in the interferogram and is usualy represented as an image where very coherent pixels have values close to one and poorly coherent pixels have values close to zero.
If a lava flow is emplaced between two successive SAR acquisitions, the interferometric coherence will be very low in the area covered by the lava flow due to the change in the geometry of the ground reflectors there. If, in contrast, the surrounding area remains coherent on the interferogram (i.e., if the lava flow is emplaced on bare rock or soil), the lava flow will appear on the interferometric coherence image as a black area surrounded by lighter-toned pixels. This is typicaly the case at Piton de la Fournaise when a lava flow is emplaced in the (largely) vegetation-free Enclos Fouqué caldera (the upper part of the volcano) where the interferometric coherence is very high (>70 %). Here, by detecting the boundary between dark and light-toned areas one can obtain a map of the lava flow contour (e.g., Rowland et al. 2003;Dietterich et al. 2012;Bato et al. 2016).
In this study, we exploited interferometric coherence images to produce an early map of the April-May 2018 lava flow. To provide relevant data for input in the response protocol in as timely fashion as possible, we used only ESA Sentinel-1 data. Sentinel-1A/B data are acquired for La Réunion island every six days both during ascending and descending passes, alternatively in Interferometric Wide Swath mode (IW, range spacing = 2.33 m, azimuth spacing = 14.06 m) and in Stripmap mode (SM, range spacing = 2.66 m, azimuth spacing = 4.15 m). The data are made freely available by the European Space Agency (ESA) via the Sentinel-1 Data Hub between 4 and 24 hours after acquisition.
As it is not possible to produce an interferogram by combining IW with SM data, so the shortest time period between two usable interferograms is 12 days. We compute the interferograms and the coherence images using the Doris 5.0.3 InSAR processor (Kampes and Usai 1999;Kampes and Hanssen 2003). To georeference the interferometric products, we used a 5 m resolution Digital Elevation Model (DEM) produced by the French Geographic Institute from two airborne LiDAR sur-veys carried out over La Réunion in 2008 and 2009. After georeferencing, the coherence maps derived from IW data have a 15 m pixel size, and those derived from SM data have a 5 m pixel size.
To discriminate, in the coherence images, areas covered by lava flow from other poorly coherent areas (e.g., due to air fall or to changes in the soil moisture between the two radar acquisitions) we developed a three step procedure. The first step consists of applying a median filter on the coherence image. Then, in a second step, we perform a clustering-based image thresholding approach: the Otsu algorithm (Otsu 1979). The resulting binary image is used to trace the lava flow boundary (including boundaries of kipukas) using the bwboundaries Matlab© function. Finally, the lava flow surface area is estimated with associated uncertainty taking into account the pixel surface and the probability that each pixel belongs to the lava flow knowing its coherence.
As the April -May 2018 eruption lasted more than one month, we were able to compute several coherence images combining, for a given acquisition geometry and mode (ascending or descending pass, IW or SM mode), a unique master image with several slave images. The master images were acquired before the beginning of the eruption and the slave images during the eruption, or just after its end. This allowed us to estimate the evolution of the lava flow surface area at different epochs (see Appendix D). Moreover, several interferograms, spanning the total duration of the eruption, were produced by combining a master image acquired before the beginning of the eruption and the slaves images acquired in the weeks following the end of the eruption. The coherence images associated with these interferograms have been stacked to produce an accurate map of the final lava flow extension and its total surface area (see Appendix D).
Generaly, the lava flow field becomes coherent a few weeks to a few months after the end of the eruption (Bato et al. 2016;Chaussard 2016;Wittmann, Sigmundsson, and Lavallée 2017). The thickness of the lava flow field can be estimated from the topographical residuals in the interferograms. These residuals reflect the deviation between a reference DEM used in the interferometric processing and the actual topography (Massonnet and Feigl 1998). Since the amplitude of topographic residuals is proportional to the perpendicular baseline of the interferometric couple (Massonnet and Feigl 1998), one can obtain a relatively accurate determination of the changes in the volcano topography due to lava flow emplacement by a statistical exploitation of several interferograms covering a large range of perpendicular baseline (Bato et al. 2016). The lava thickness obtained in this way can be used not only to evaluate the volume of the lava flow field, and then, by making some assumptions on the lava porosity, the volume of lava emitted, but also to update the DEM after each eruption. Such an update is mandatory to achieve model-based flow-length projections. Unfortunately, in the case of the April 2018 eruption, we could not use the Sentinel-1 interferograms to determine the lava flow thickness since the perpendicular baseline (B perp ) of the S1 interferometric couples was always very low ( 48 m 35 m). A particular effort has been made by ESA in the Sentinel-1 system design to achieve an orbital tube of 50 m radius (rms). This guaranties high performances for ground surface displacement measurement purposes by reducing the sensibility of the interferograms to possible topographic artifacts. However, in our case, this "improvement" is a disadvantage since we are looking for topographic residuals (Ebmeier et al. 2012;Albino et al. 2015;Kubanek, Westerhaus, et al. 2015;Kubanek, Richardson, et al. 2015;Bato et al. 2016;Kubanek, Westerhaus, and Heck 2017). The characteristics of CosmoSky-Med and TerraSAR-X/TanDEM-X constellations make them better suited for this type of application. The use of future CSK and TSX/TDX acquisitions for Piton de la Fournaise will allow us to calculate the April-May 2018 lava flow thickness. atellite S1B S1B Length = S1A S1A Length = S1B S1B Length > S1A S1A Length = d flow area m atellite S1B S1B