Monitoring the time-averaged discharge rates, volumes and emplacement style of large lava flows by using MIROVA system: the case of the 2014-2015 eruption at Holuhraun (Iceland)

25 The 2014-2015 eruption at Holuhraun has produced more than 1.5 km 3 of basaltic magma and can be 26 considered one of the major effusive events of the last two centuries in the world. During this eruption 27 the MIROVA system (a volcanic hot-spot detection system based on MODIS middle infrared – MIR 28 - data) has been used to detect and locate the active portions of the lava flow(s), and to measure the 29 heat radiated by the growing lava field. According to these data the eruption was characterized by a 30 slow decay of the radiant power, accompanied by a change in the lava transport mechanism that 31 shifted from open channels, at the beginning of the eruption, to lava tubes, during the last months of 32 activity. Despite the evident evolution of lava transport mechanism, we found that the overall 33 decreasing trend of the thermal flux was mainly controlled by the exponential decline of lava 34 discharge rates, while the increasing insulation of the flow field had a strong impact in transporting 35 efficiently the lava at the distal flow front(s). Our results suggest the apparent time averaged lava 36 discharge rates ( TADR ), derived from satellite thermal data, may fluctuate around the real effusion 37 rate at the vent, especially in the case of large lava flows emplacing in nearly flat conditions. The 38 magnitude and frequency of these fluctuations are mainly controlled by the emplacement dynamic, 39 (i.e. occurrence of distinct major flow units), while the transition from channel- to tube-fed lava 40 transport mechanism play only a minor role (±30%) in the retrieval of TADR using MIR data . When 41 the TADR values are integrated to calculate erupted lava volumes, the effects of pulsating 42 emplacement dynamic become smoothed and the eruptive trend become more clear. 43 We suggest that during the Holuhraun ’s eruption, as well as during many other effusive eruptions, 44 the MIR-derived radiant flux essentially mimic the overall trend of lava discharge rates, with only a 45 minor influence due to the emplacement style and evolving eruptive conditions. From a monitoring 46 and operational perspective, MIROVA demonstrates to be a very valuable tool to follow and, 47 possibly, forecast the evolution of on-going effusive eruptions.


INTRODUCTION
On August 29 th 2014, one of the largest effusive erup− tion of the last 3 centuries took place along the Eastern Volcanic Zone (EVZ) of Iceland (Figure 1), about 45 km north−east of Bárðarbunga volcanic system [Gudmun− sson et al., 2014].The new eruption followed 15 days of sustained seismicity and accompanied the propagation of a 45 km−long dike [Sigmundsson et al., 2015] that un− locked an historical eruptive path named Holuhraun [Sig− urðsson and Sparks, 1978].The emission of lava persisted at very high rate for 180 days up to February 27 th , 2015 when the activity was declared over [Gislason et al., 2015].The eruption was characterized by a slow decay of ef− fusion rate [Coppola et al., 2017] that accompanied a co− eval slow subsidence of the Bárðarbunga caldera [Gud− mundsson et al., 2016].The clear link between the two pro− cesses provided clear evidences of an inelastic eruption, whereby the lateral magma withdrawal was essentially linked to the gravity−driven collapse of the summit caldera [Coppola et al., 2017;Gudmundsson et al. 2016].About 84 km 2 of nearly flat land was covered by the new lava field that reached a maximum extent of about 18 x 5 km (Figure 1) and thickness up to 40 m [Gislason et al., 2015;Pedersen et al. 2017].
The prohibitive environmental conditions that might characterize these latitudes made field observations at the eruptive site very hard, especially on a continuous basis and for several months during the winter time.In these conditions, space−based thermal data have been extremely useful since they allowed a safe detection, location and quantification of the radiant flux produced by the effu− sive activity.In particular, infrared data acquired by MODIS (Moderate Resolution Imaging Spectroradiometer) and elab− orated through the MIROVA volcanic hotspot detection sys− tem, [Coppola et al., 2016;www.mirovaweb.it],were de− livered to the Icelandic Meteorological Office (IMO) as part of its daily operational monitoring activity of the Holuhraun eruption.
Since the beginning of the eruption the main tasks of satellite thermal data were: (i) to provide information about the location of active lava flow areas and front, (ii) to give an estimation of lava discharge rates and erupted volumes and finally (iii) to identify the ongoing effusive trend (steady, waning or waxing).The latest two points became a priority after a couple of months when the winter sea− son and the low−light conditions made the eruptive site hardly accessible and field observations became more and more rare and difficult.The calculation of effusion rates from satellite thermal data (the so called "thermal proxy"), become one of the main topic discussion among a ther− mal remote sensing group, established in that occasion by Dr. Thor Thordarson [Harris et al., 2016a].One of the main questions discussed within the group was "whether the sub− tle, but persistent decay of thermal flux indicated by satel− lite measurements, was caused by a real, slow decline of the effusion rate or rather it reflected the increasing in− sulation of the growing lava field and the development of a lava tube system [Pedersen et al., 2017].
In this paper we present time−averaged lava discharge rates and volumes, derived by using MIROVA system, and we compare the results with field and independent mea− surements collected during and after the Holuhraun erup− tion.The comparison reveals the role of changing em− placement style in the lava discharge rates calculation and outlines the contribution of MIROVA in safely tracking this large effusive eruptions from space.

TIME-AVERAGED DISCHARGE RATES FROM SATELLITE THERMAL DATA: BACKGROUND AND OPEN QUESTIONS
The relationship between effusion rates and thermal emissions of lava flows has received increasingly at− tention during the past three decades (a full list of 46 papers published between 1990 and 2005 is reported in Harris [2013]).The methods, limits and applications of this approach are part of a book expressly focussed on detecting, modelling and responding to effusive erup− tions [Harris et al., 2016b], whereas exhaustive overview of the physic behind mass and energy flow through a lava flow system is provided by a series topical works [Pieri and Baloga 1986;Crisp and Baloga 1990;Harris et al., 1997;Wright et al., 2001;Harris et al., 2007;Har− ris and Baloga, 2009;Dragoni and Tallarico, 2009;Cop− pola et al., 2013;Garel et al., 2012Garel et al., , 2014Garel et al., , 2015;;Tarquini 2017;among others].Here, we outline the basic princi− ples of this technique and we summarize the open ques− tions that have been addressed in this work.
To describe volumetric flow rates of erupted lavas, we use the terminology given by [Harris et al., 2007].We use the generic term effusion rate, to describe the instanta− neous rate at which lava is erupted from the vent at any time.The term mean output rate (MOR) is used to de− scribe the final volume of erupted lava divided by the total duration of the eruption.Finally, we used the term time−averaged lava discharge rate (TADR) to describe the volume of lava emplaced during a specific time in− terval [Harris et al., 2007].This definition better applies to the use of satellite−based methods, for measuring the changes in lava volume over a given period of time prior the image acquisition [Wright et al., 2001].
Currently, two main methods exist to estimate TADRs from space−based thermal data; the method of Harris et al., [1997], simplified later by Wright et al. [2001], and the method of Coppola et al., [2013].Both the methods rely on the original heat balance approach [Pieri and Baloga, 1986], stating that the mean output rate MOR (m 3 s −1 ) of a cooling−limited lava flow, is related to its final plan area A (attained when the flow achieves its fi− nal length, L), by: (1) where ρ (kg m −3 ) and C p (J kg −1 K −1 ) are the bulk den− sity and heat specific capacity of lava, respectively, T e (K) is the effective radiation temperature of the flow, T 0 (K) is the lava eruption temperature, and T stop (K) is the temperature of the flow front at the time the flow has cooled to a halt.In this framework, the effective radi− ation temperature is defined as "the temperature at which the flow would radiate if it had a constant sur− face temperature throughout the emplacement of the flow" [Crisp and Baloga, 1990].Harris et al., [1997], applied this approach to satel− lite thermal data, developing what is actually called the "thermal proxy".Further works [Wright et al., 2001;Harris et al. 2007Harris et al. , 2009] ] refined and simplified this method suggesting that the equation 1, when applied to satellite radiance data, provides discharge rates that are not necessarily averaged over the total duration of an eruption, but rather over "some time" prior to the satel− lite acquisition.Here, the effective radiation temperature, T e , takes a specific meaning, since it represents the av− erage surface temperature over "some time" prior to the satellite acquisition, weighted according to the radiative heat flux [Harris and Baloga, 2009].Accordingly, Har− ris et al. [2007] proposed to use the term "time−averaged discharge rate" (TADR), in order to describe the discrete measurements of lava flux retrieved from single satel− lite thermal images.Wright et al. [2001] also argued that all values except MOR (or TADR) and A are assumed a priori in equation 1 so that the thermal proxy reduces to a simple relationship whereby: (2) Here x is an empirical parameter (m s −1 ), that de− scribes the appropriate compositional flow parameters (ρ, C p ) as well as thermal insulation (T e ) and cooling (T 0 −T stop ) conditions [cf.Harris and Baloga, 2009].In this formulation the area of the active flow area, A (m 2 ) is also dependent on the insulation condition expressed by T e , and is directly retrieved from the pixel integrated spectral radiance, R λ (W m −2 sr −1 μm −1 ), according to: (3) where T bk is the background temperature (K), L λ is the Plank function for wavelength λ (W m −2 sr−1 μm −1 ), and Apixel is pixel area (1 km 2 for MODIS).By assum− ing two end−member radiating temperatures (T e ) for the hot and cold models (for example 500°C and 100°C for the channel− and tube−fed flow type, respectively), the equation 3 allows the user to calculate a range of plau− sible active flow areas, responsible for the observed ra− diance [Harris et al., 2007].
As stressed by Harris and Baloga [2009], the values for x have to be set on a case−by−case, thus leaving a wide arbitrariness to the user that may adjust all the un− knowns in equation (1) to achieve a best-fit with avail− able and independent field data.The TIR (thermal in− frared) bands are generally used in this approach under the assumption that the surface area of high tempera− ture cracks is small and provides a negligible contribu− tion to the pixel integrated radiance at 10 μm to 12 μm [Harris and Baloga, 2009] This approach (equations 2 and 3) was the one used by Bonny et al. [2018] to esti− mate the volume of lava erupted during the 2014−2015 eruption at Holuhraun by using MODIS TIR data.
The method of Coppola et al. [2013], does not hold with the calculation of active flow areas, but is based on the empirical relationship that directly links the TADR with the Volcanic Radiant Power (VRP) calculated via the MIR method [Wooster et al., 2003].The MIR−derived VRP is a measurement of the heat flux radiated almost exclusively by the portions of lava flows having effec− tive temperature (T e ) above 600 K, and is calculated as: (4) where R MIR,alert is the pixel integrated MIR radiance of the i th alerted pixel, R MIR,bk is the MIR radiance of the background, A pixel is the pixel size (1 km 2 for the re− sampled MODIS pixels), and 18.9 is a constant of pro− portionality [Wooster et al., 2003].The TADR is than calculated by using a single coefficient called radiant density (c rad in in J m −3 ): (5) that describes the relationship between volumetric and radiant flux appropriate for the observed eruption [Cop− pola et al., 2013].By analyzing a large compositional spectrum of recent lava flows, Coppola et al. [2013] pro− posed an empirical method to calculate the c rad , based on the silica content of the erupted lavas: (6) where X SiO2 is the silica content of the erupted lavas (wt%).
This method allows to take into account the strong effects that the bulk rheology has on the spreading and cooling processes of active lavas, and, according to the authors, permit to estimates TADR with an uncertainty of ± 50% [Coppola et al., 2013].Most importantly, it al− lows to calculate the appropriate radiant density, once the first chemical analysis of the erupted lavas become available [Coppola et al., 2017].Since the MIROVA sys− tem provides automatically the VRP [Coppola et al., 2016], the radiant density approach (equations 5 and 6) was the one used by Coppola et al. [2017] to estimate TADR and erupted lava volume during the Holuhraun eruption.
Although the two methods derive from the original heat balance approach, valid for cooling−limited lava flows (Equation 1), it is important to remark that they do not necessarily yield exactly the same results [cf.Bonny et al., 2018;Coppola et al., 2017].This is because the wavelengths used by the two methods are different (R TIR and R MIR , respectively) but also because of the arbitrari− ness of the input coefficients, as x for the Harris' method and c rad for the Coppola's method.Nevertheless, both the methods provide an upper and lower boundary estimates of TADR that are linearly related to the pixel−integrated spectral radiances (cf.Equation 1 to 5).It follows that whatever the method used, a constant TADR would be translated into a constant radiance detected by space.
This point was challenged by recent laboratory ex− periments and physical modelling of cooling viscous gravity currents [Garel et al., 2012[Garel et al., , 2014[Garel et al., , 2015] ] accord− ing to which for a given magma discharge rate the heat radiated by the flow surface reaches a steady value only after a transient time.According to the cooling− limited flow model (Equation 1), this transient time would correspond to the time required for the lava flow front to reach its maximum length (L) and to cool to a halt (T front = T stop ).Hence, during rapid changes of ef− fusion rates the thermal signal would be "buffered" in time, because it takes time for the higher flow rate to propagate downstream and cause perceptible increases in flow area [Harris and Baloga, 2009].The transient time thus reflect the dynamic response of the thermal structure toward the appropriate cooling−limited con− ditions expressed by the new TADR vs.A relationship (i.e.new steady−state condition).Tarquini [2017] also suggests that this transient time is typical of non−equi− librium steady−state systems, and is related to a struc− tural relaxation time that operate particularly on long− lived lava flows.In this non−equilibrium framework, the active lava flow units represent "dissipative structures" that promote significant fluctuations in the radiative power, even in the case of a constant supply.Hence, the radiant flux would be modulated by a pulsating em− placement dynamic, associated to the occurrence of flow diversions that resets the system further from its maximum extension [cf.Tarquini, 2017].
Other open questions remain challenging, especially for operational use of the thermal proxy during effusive crisis.For example, based on its theoretical model, Garel MONITORING DISCHARGE RATES OF LARGE LAVA FLOWS et al. [2012] suggests that when lava tubes form, the sur− face thermal signal will not reflect the flow dynamics, be− cause the low crust temperatures do not reflect a poten− tially high flow rate of hot lava underneath.According to the authors, this would pose a serious problem, since many lava flows exhibit evolving emplacement styles throughout the eruption, often characterized by increas− ing insulation conditions and gradual formation of lava tubes.Similarly, it is still unclear how to determine whether an eruption has ceased based on thermal data, or whether the lava stored underneath a cooled crust can potentially emerge suddenly and spread rapidly [Garel et al., 2015].In these terms the 2014−2015 eruption at Holuhraun represents a key−case since allowed us to address all these points over the largest effusive eruptions occurred in the past three centuries.

THE MIROVA SYSTEM
MIROVA (Middle Infrared Observation of Volcanic Activity) is an automated global hot spot detection sys− tem (www.mirovaweb.it)based on near−real time pro− cessing of Moderate Resolution Imaging Spectrora− diometer (MODIS) infrared data [Coppola et al., 2016].The system completes automatic detection and location of thermal anomalies, and provides a quantification of the Volcanic Radiative Power (VRP), within 1 to 4 hours of each satellite overpass.This is achieved through a hy− brid algorithm [fully described in Coppola et al., 2016] based on MIR radiance data (~ 4 μm) recorded at 1−km spatial resolution.Two MODIS instruments (carried on two NASA's satellites: Terra and Aqua) deliver approx− imately 4 images per day for every target volcano lo− cated at equatorial latitudes.However, due to the polar, sun−synchronous orbit, the two satellites increase their sampling time at high latitudes, providing, at least 6 to 10 overpasses per day over Iceland.An example of thermal images elaborated by the MIROVA system dur− ing the Holuhraun eruption is given in Figure 1b.

VOLCANIC RADIATIVE POWER AND ENERGY OF THE 2014-2015 ERUPTION AT HOLUHRAUN
Between August 29 th , 2014 and March 4 th , 2015, MIRO− VA detected hotspots in 1105 images over a total of 1623 MODIS overpasses above Iceland (~68.1%).The volcanic radiative power (VRP) ranged from ~39.1 GW, during the initial stage of the effusion, to less than 10 MW just be− fore the end of the eruption (Figure 2a).During the course of the eruption, we visually inspected all the acquired scenes and we identified a large number of images acquired in cloudy conditions, and/or under poor geometrical con− ditions that strongly deformed and affected the thermal anomaly at ground level.These images were discarded from further analysis (blue circles in Figure 2a) so that a su− pervised dataset of only 206 high−quality images (~12.7% of the total MODIS overpasses) was used to pro− vide a robust quantification of VRP produced by the erup− tion (red circles in Figure 2a).As a whole we estimated that the Holuhraun eruption radiated approximately 1.6 ×10 17 J into the atmosphere (red curve in Figure 2b), with a time− averaged radiant flux of about 10.5 GW.

TIME-AVERAGED DISCHARGE RATES AND VOLUME CALCULATION
During the course of the eruption, we calculated TADRs and erupted lava volumes by using the radiant density approach described above (Equation 5 and 6).The pre− liminary chemical analysis of Holuhraun lavas (SiO 2 = 50.5 wt%), provided in early September 2014 by the Uni− versity of Iceland (http://earthice.hi.is/bardarbun− ga_2014), were used to calculate (Equation 6) a radiant density comprised between 0.62 × 10 8 and 1.86 × 10 8 J m −3 .These two values allowed us to provide an upper and lower boundary limits for lava volume calculation, with the mean value being the most likely (Figure 3a).Given a final VRE equal to 1.64×10 17 J (Figure 2b), this method provided a first, timely assessment of the erupted lava volume equal to 1.75 ± 0.88 km 3 (Figure 3a).This was in good agreement with the flow volume es− timated a−posteriori ( 1.4 − 1.6 km 3 ) by using indepen− dent datasets [Gislason et al., 2015;Gudmundsson et al., 2016;Dirshel and Rossi, 2018;Bonny et al., 2018].The inferred satellite−derived TADR and erupted volumes were updated continuously and delivered via emails to the teams of the IMO and University of Iceland in charge of the monitoring.

ERUPTIVE TREND DERIVED FROM SATEL-LITE THERMAL DATA
One of the most significant feature of this eruption was the overall exponential trend of the satellite−derived cumulative volume (Figure 3b).As outlined by recent studies, this trend mirrors almost perfectly the slow subsidence of the Bárðarbunga caldera floor (Figure 3b) and suggests a strong connection between source caldera/reservoir dynamic, and the eruptive dynamics 45 km distant [Gudmundsson et al., 2016;Coppola et al., 2017].According to these works, the exponential decay of effusion rates was driven by a decrease of magma− static pressure, associated to the shrinking of the magma chamber and the gravity−driven collapse of Bárðar− bunga caldera.In this scenario, the instantaneous effu− sion rate at the vent follows a perfectly exponential curve (blue line in Figure 4a) that has been used as a benchmark for our MODIS−derived TADRs' (red line in Figure 4a).Hence, the residuals between the two rates (observed − modeled; Figure 4b) enhance temporary variations of the source process, or instabilities of the thermal proxy due to lava flow emplacement dynamic.In the following section, we discuss the effusive trend in the light of the evolving emplacement style observed in the Holuhraun lava field.
According to Pedersen et al. [2017] the Holuhraun eruption can be subdivided in three main phases, char− COPPOLA ET AL. acterized by evolving lava transport processes.
The phase 1 (31 August − 12 October 2014) was char− acterized by lava transport confined to open channels that led the formation of four consecutive lava flows emplaced side by side (no.1−4 in Figure 4b).During this phase the discharge rates measured on the field were comprised be− tween 350 and 100 m 3 s −1 [Pedersen et al., 2017], in good agreement with satellite−derived values.The four flow units reached distances of 16.9, 11.8, 16.3 and 11.0 km, respectively (Figure 4a), and covered a total area of 58.31 km 2 [Pedersen et al., 2017].Given a flow field volume of 0.78 km 3 estimated on October 12 (Figure 3b), the mean flow thickness characterizing the Phase 1 was 13.4 m, consistent with the thickness of each single flow units ob− served on the field [Pedersen et al., 2017].During this pe− riod our TADRs measurements show the greatest oscilla− tions around the modelled value (up to ± 60%), resulting in sharp variations exactly in correspondence of the ac− tivations of the distinct flow units (Figure 4b).
The phase 2 (13 October − 30 November 2014) had lower effusion rates (150 − 100 m 3 s −1 ) and was charac− terized by the formation of a < 1km 2 lava pond, that acted as a distributor of subsequent lava flows [i.e.flows no.5− 8; Pedersen et al., 2017].These flow units reached grad− ually shorter distances (11.7, 8.4, 7.1, 5.4 km, respectively) and were accompanied by satellite−derived TADRs show− ing decreasing oscillation amplitudes (Figure 4b).Towards the end of this phase the occurrence of inflation plateaus provided the first evidences of a growing lava tube sys− tem, with satellite−derived TADRs becoming more stable (Figure 5b).About 0.52 km 3 , of lava was erupted during this phase (Figure 3a) producing a net increment of the flow field area of 18.7 km 2 [Pedersen et al., 2017].By the end of November the flow field reached an average thick− ness of 16.9 m.
During the third phase (1 December 2014 − 27 Febru− ary 2015) the lava transport was mainly confined to lava tubes which fed several breakouts and inflation plateaus along the initial flow units no. 1 and 2 [Pedersen et al., 2017].More than 19 km 2 of the flow field was resurfaced, with discharge rates much more steady (±20%) and typ− ically lower than 100 m 3 s −1 .A sharp reduction of TADR was recorded since 27 January 2015 and preluded the end of the eruption occurred one month later, on the 27 February 2015 (Figure 4).This reduction was interpreted by Coppola et al., [2017] as due to the gradual closure of the magma path once the overpressure inside the dike had dropped below a critical value.The lava tube system distributed about 0.4 km 3 during this phase, with the lava field reaching a final extension of 85.4 km 2 .The infla− tion and resurfacing processes that operated through the tube system led the average flow thickness to increase up to 20.2 m.
The emplacement of the major lava flows (phases 1 and 2) is outlined by the occurrence of apparent TADR's "pulses" (Figure 4a) oscillating around the modelled value for effusion rate (Figure 4b).Pulsating effusive activity has been recognized at several basaltic volcanoes and may results from different processes such as pulsed magma supply, repeated accumulation and collapse of a foam layer at the reservoir roof, or by processes of magma mixing occurring within the magma chamber [e.g., Har− ris and Neri, 2002;Lautze et al., 2004].However, in the case of Holuhraun eruption, the link between the effusive trend and the collapse of the caldera (Figure 3b) claims for a simple, gravity−driven dynamic [Gudmundsson et al., 2016;Coppola et al., 2017] that would exclude, or minimize, the occurrence of above mentioned processes.Hence, the TADR's pattern overprinted to the exponen− tial trend is likely related to the emplacement dynamic occurred during the course of the eruption, and especially during the phases 1 and 2. Laboratory experiments [Garel et al., 2012[Garel et al., , 2014] ] and theoretical treatments [Tarquini, 2017], demonstrated that a constant effusion rate can ac− tually produce pulses of the radiative power (Figure 4b) that reflect the non−equilibrium steady−state growth of compound lava flow fields.According to Tarquini [2017] each flow unit will advance until the thermal equilibrium will be reached and the flow front stop.Hence, the sys− tem is resettled to give origin to a new flow unit.This seems to be exactly the dynamic occurred during the em− placement of the Holuhraun lava field, in which each ap− parent TADR's pulse correlates a phase of lengthen− ing of a distinct flow unit.Notably, the amplitude of the pulses decreased with time as the overall effusion rate de− clined, the maximum length attained by the flow units re− duced, and the emplacement style evolved from channel− to tube−fed (Figure 4b).

EVOLVING EMPLACEMENT STYLE, EVOLVING THERMAL STRUCTURE
In this section, we illustrate how the evolving em− placement style and the channel− to tube−fed transition affected the surface thermal structure of the Holuhraun's lava field.In particular, we used the temperature profiles (as the one shown in Figure 1c) obtained from the MODIS image, in order to track the position of the active flow fronts and all the hot lava surfaces emerging along the principal flow path (from West to East).By stacking to− gether all the W−E profiles the presence of active lavas (i.e. the vent(s), open−channel(s), breakout(s), flow front(s), etc.) appears as high−temperature features (pix− els), whose intensity, contiguity and extension can be easily visualized as shown in Figure 5a.We assumed that the position of the flow front is represented by the east− ernmost pixel having a brightness temperature (BT4) of 50°C higher than the background (white dashed line in Figure 5a).Hence, the distance between this pixel and the one located over the active vent (Figure 1b), provides a measurement of the active flow length at each satellite overpass (red line in Figure 5b).Lava flow lengths mea− sured on the field during the phases 1 and 2 (Flows no. 1 to 8) are also plotted for comparison (Figure 5b).This analysis indicates that during the phases 1 and 2, the active flow fronts were thermally connected to the vents, through a high−radiating feature representing the open channels.On the other hand, since December 2014 the reduction of lava discharge rate caused a marked shortening of the open channels, although the lava field was still characterized by the presence of ac− tive lavas at the distal front(s).Notably, between the two hot, active zones (i.e.open channels and flow front) there is a low−temperature zone corresponding to crusted and cooling lava surfaces (Figure 5a).This is a clear thermal evidence that since the beginning of the phase 3 the transport of lava occurred along a master tube system and that the emplacement style started to be tube−fed dominated.
As a pure theoretical exercise, the flow extension measured by using the temperature profiles has been compared (Figure 5b) with the maximum distance (L max ) predicted by one of the simplest empirical model de− veloped for single, cooling−limited, channel−fed lava flow [cf.Walker, 1973;Harris and Rowland, 2009].In particular, we used the Calvari and Pinkerton's equation [1998] that takes into account exclusively the TADRs measurements for the estimation of flow length: (7) The coefficients of equation 7 were empirically de− termined by regression analysis of single channel−fed Etnean lava flows lasting 1 to 12 days, a duration very similar to the single flow units observed at Holuhraun.This model does not require any knowledge of several others parameters (i.e.underlying slope, lava tempera− ture, viscosity, velocity, channel width, channel depth, etc…) that are necessary to predict flow length using other approaches [see Harris and Rowland, 2009 for a review], but that are difficult, if not impossible to be constrained for each flow unit.However, despite its simplicity, this empirical model of Calvari and Pinker− ton 1998 allows to calculate the maximum length of the Holuhraun's lava flow, and to see the effects of an em− placement style persistently dominated by channel−fed, cooling−limited conditions (blue line in Figure 5b).
During the phases 1 and 2, there is an excellent cor− relation between the measured and the modelled flow lengths, both in terms of absolute values and pattern (Figure 5b).This corroborate the fact that during this pe− riod the emplacement style was effectively those for which equation 7 was proposed to describe, that is sin− gle, channel−fed flows.Differently, since December 2014, there is a growing discrepancy between the mea− sured and modelled flow length (Figure 5b), exactly in correspondence of the formation of a master lava tube, as observed in the field [Pedersen et al., 2017].
From this picture, it is clear that the lower discharge rates characterizing the Phase 3, did not affect the max− imum extent of the lava flow, still reaching a distance of more than 15 km, just before the end of the eruption (Figure 5b).By paraphrasing Harris and Rowland [2009] "the transition from poorly insulated (channel−fed) to well insulated (tube−fed) regimes increases the length that a lava flow can extend for a given TADR".
Notably, we may calculate the radiant density (c rad = VRE/Vol) typical of the tube−fed flow field, by di− viding the thermal energy measured during the Phase 3 (VRE = ~ 0.4 × 10 17 J; Figure 2b) by the volume lava erupt− ed in the same time window, derived from Tandem−X mea− surements (Vol = 2.6 × 10 8 m 3 ; Figure 3a) [Dirscherl and Rossi, 2018].The resulting radiant density (c rad = ~1.5 × 10 8 J m −3 ) is only slightly higher than the radiant den− sity calculated in the same way for the bulk channel−fed flow field emplaced during phases 1 and 2 (c rad = ~1.1 × 10 8 J m −3 ).This suggests that the variation of the em− placement style produces an uncertainty in the TADRs measurement that does not exceed 30%, when provided in near real time.In other words, our results point out that, the formation of a lava tube does not affect substantial− ly the amount of "active" lava radiating to the atmosphere, but simply "made this lava to be transported further from the vent".The lava tube thus act as a simple extension of the volcanic conduit, by displacing the vent down− stream, likely along the main axis of the previously formed flow units.The recognition of this behavior in near real time allow to update the vent location and has been proved to be a necessary improvement to correct estimation of the lava flow path and the effective runout distance [Harris et al., this issue].

OPERATIONAL USE OF MIROVA SYSTEM DURING THE EFFUSIVE CRISIS
By the end of November 2014, MIROVA started to be part of the monitoring system operative at the Ice− landic Meteorological Office (IMO) with the coverage of nine active volcanic systems.IMO is the main institu−

MONITORING DISCHARGE RATES OF LARGE LAVA FLOWS
tion in Iceland in charge of monitoring natural hazards and issuing forecast and warnings in case of impend− ing or on−going eruptions.It operates a multi−sensors network that includes seismometers, GPS and gas sen− sors for the geophysical monitoring.During the Holuhraun eruption, gas monitoring activity was im− proved and new instrumentations were deployed and maintained at the eruption site for the entire duration of the event [Pfeffer et al., 2017].The MIROVA system showed its strength and reliability in a period of the year when, in Iceland, the day light time is very short.Sev− eral of the instrumentations and products that could normally be used to monitor an eruption, started to fail or to give incomplete data sets due to harsh weather con− ditions and limited sunlight.MIROVA provided a robust and continuous data set on a daily basis which was used from November to March to follow the health and strength of the eruption, as well as the post−eruption phase.
On the operative side this system is revealed to be an important monitoring tool for following the temporal evolution of the eruption and providing estimates of the erupted lava volume in near real time [Coppola et al., 2017].Extrapolation of the effusive trend during the on− going eruption may in fact be used to forecast total erupted lava volume, which in turn is fundamental to drive some lava flow modelling [i.e.Tarquini et al., this issue].Also, the recognition of emplacement style through the analysis of flow' surface thermal structure (i.e.channel− or tube−fed lava flows) can be used to up− date the position of active flow fronts, and to choose the appropriate modeling strategy, by tuning some appro− priate parameters that govern lava flow simulation's codes [i.e.Tarquini et al., 2018, this issue].
Most importantly the thermal data have proved to be very important for constraining the time of the eruption end.A drastic change of TADR was observed since the end of January 2015 becoming clear the second week of February 2015 (Figure 4a).This change was interpreted by Coppola et al. [2017] as due to the gradual closure of the magma path (dike) and indicated the potential tim− ing of the end of the eruption in quite a good advance, about one month before the declaration of ceased erup− tion.These data were discussed within the Scientific Board Management meetings that were held from the beginning of the eruption in collaboration with the Icelandic Civil Protection and the University of Iceland.The very low TADR detected on 26 February 2017 (0.2 m 3 s −1 ; Figure 4a), supported the decision for a field crew to flight over the eruptive site on February 27 and de− clare that the eruption was over.
The thermal data were also useful for providing a semi−quantitative indication of the strength of the gas source affecting the ground level concentration [Sim− mons et al., 2017].Due to its extension the lava field it− self was a considerable source of volcanic gases.In ad− dition, the gas released close to the surface was more easily trapped within the boundary layer and largely af− fected the concentration at ground.This fact suggested that the MIROVA detection and the VRP temporal vari− ation could be correlated with the amount of gases re− leased at the source.In this way the VRP data have been used by the forecasters at IMO to identify, in a qualita− tive way, when to expect an increase in the SO 2 values in areas located downwind the eruption site.

CONCLUSIONS
During the 2014−2015 Bárðarbunga−Holuhraun eruptive crisis, one of the main challenge for the vol− canological community was to monitor the effusive process in a safely, timely and routinely manner.For the whole duration of the eruption, MIROVA provided a ro− bust and continuous set of thermal data that were used to depict the effusive trend and to estimates the erupted lava volume although the changing emplacement con− ditions (from channel− to tube−fed flow) introduced un− certainties in the interpretation of the data.
The retrospective analysis of this eruption provides an exceptional opportunity to study the relationship be− tween heat flux and TADR during the emplacement of a large compound lava field.We show that the overall trend of thermal emission was related to the combined effect of two over imposed patterns: (i) a main expo− nential decay of the effusion rate trend, governed by the source processes and (ii) a secondary pulsating pattern related to the emplacement dynamic of the flow field.We found that the magnitude and timing of these pulses were strictly related to the timing and length scales of discrete flow units characterizing open−channel lava transport.Conversely, the formation of lava tubes pro− duces smaller instabilities and promotes lava to flow at greater distance in steady−state thermal conditions.This process is clearly visible from the evolution of the thermal structure depicted from the MODIS−MIROVA images and allow to track changes in the lava transport mechanisms operating during the eruption (Figure 5a).

COPPOLA ET AL.
The results presented here clearly indicate that the cal− culation of erupted lava volumes from the integration of satellite−derived TADR allow to smooth the short−term perturbation associated to the emplacement dynamic and provide a robust way to depict source eruptive trends.We thus regard the application of this methodol− ogy as a key−factor in volcano monitoring and satel− lite-data-driven response to an effusive crisis [Harris et al., 2018].
MIROVA is a collaborative project between the Universities of Turin and Florence (Italy) and is supported by the Centre for Volcanic Risk of the Italian Civil Protection Department.We acknowledge the LANCE-MODIS data system for providing MODIS Near Real Time products.

FIGURE 1 .
FIGURE 1.(a) Location of the Bárðarbunga−Holuhraun volcanic system showing the ice covered, 10−km−wide caldera of Bárðar− bunga and the 2014−2015 lava filed at Holuhraun.(b) Example of MODIS−derived thermal image (spatial resolution of 1 km) showing the Brightness Temperature at 4 μm (BT4) recorded on February 03rd, 2015 over the Holuhraun lava field (white line).The star indicates the approximate location of the vent.(c) Temperature profile along W−E direction, ob− tained by calculating the maximum BT4 along each vertical (N−S) line.The position of the flow front is marked by the easternmost thermally anomalous pixel (see the text for details).

FIGURE 2 .
FIGURE 2. (a) Timeseries of VRP retrieved from MODIS−MIROVA system during the 2014−2015 eruption at Holuhraun.Blue cir− cles represent data discarded because of cloudy conditions or poor viewing geometry.The data selected for VPR calcu− lation are the white circles.The red line is a 5−points−moving average of VRP.A VRP below 10 MW coincides with the end of the eruption.(b) Volcanic Radiative Energy (VRE) calculated by trapezoidal integration of selected data only (red curve) or all the data (blue curve).Note the strong influence of data selection on the final total energy.Approximately 1.6×10 17 J were radiated into the atmosphere (red curve), with a mean radiant flux of about 10.5 GW.

FIGURE 3 .
FIGURE 3. (a) Lava volumes measured during the 2014−2015 eruption at Holuhraun.The colored field envelop the upper and lower boundary limits calculated using the radiant density approach (eqs.5 and 6).Best−estimate volume calculation (blue cir− cles), and relative uncertainty (±50%) were provided in near real time to IMO and University of Iceland.Independent measurements of final lava flow volumes given by several authors are also shown.(b) Comparison between lava field volume measured by MODIS (blue circles) and exponential models of erupted lava volume (Coppola et al., 2017) and caldera volume [Gudmunsson et al., 2016].The good correlation indicates a strong link between to the gravity−driven collapse of Bárðarbunga caldera and the gradual decrease of magma−static pressure driving the effusive eruption at the vent [cf.Gudmunsson et al., 2016; Coppola et al., 2017; Dirsherl and Rossi, 2018].

FIGURE 4 .
FIGURE 4. (a) derived from MIROVA data (white circles) by using the radiant density approach (eqs.5 and 6).The blue line represents the best−fit exponential model [Coppola et al., 2017] hereby considered the real effusion rate at the vent.(b) Residuals flow rates (left axis) obtained by subtracting the model from the observations (i.e.MODIS−derived TADR minus mod− elled effusion rates).Note how the residuals oscillate around the model showing a decreasing amplitude through time.The timing and amplitude of the resid− ual pulses are coherent with the occurrence and maximum length of the main channel−fed flow units numbered from 1 to 8 [data from Pedersen et al., 2017].The three phase of the eruption described by Pedersen et al. [2017] are also shown.

FIGURE 5 .
FIGURE 5. (a) Evolution of the along−path thermal structure of the Holuhraun lava flow.This plot is obtained by stacking together all the thermal profiles (max BT4 along each N−S line) as those illustrated in Figure 1c.The white dashed line indicate the position of the active flow front (see the text for details); (b) Comparison between measured (red dashed line) and modelled (blue line) flow length.The Calvari & Pinkerton's model (eq.7) assumes emplacement style dominated by chan− nel−fed, cooling limited conditions.The discrepancy between measured and modelled flow distance after December 2014 outline the increasing efficiency of tube−fed emplacement style in transporting the lava to the flow front.