“MODELING LAVA FLOW PROPAGATION OVER A FLAT LANDSCAPE BY USING MRLAVALOBA: THE CASE OF THE 2014–2015 ERUPTION

During the emplacement of the 2014−2015 lava flow in Holuhraun (Iceland) a new code for the simulation of lava flows (MrLavaLoba) was developed and tested. MrLavaLoba is a probabilistic code which derives the area likely to be inundated and the thickness of the final lava deposit. The flow field in Holuhraun progressed through a fairly flat floodplain, and the initial limited availability of topographic data was challenging, forcing us to develop specific modeling strategies. The development of the code, as well as simulation tests, continued after the end of the eruption, and latest results largely benefitted from the availability of improved topographic data. MrLavaLoba simu− lations of the Holuhraun scenario are compared with detailed observational analyses derived from the literature. The obtained results high− light that small−scale morphological features in the pre−emplacement topography can strongly influence the propagation of the flow. The distribution of the volume settling throughout the extension of the flow field turned out to be very important, and strongly affects the fit between the simulated and the real extent of the flow field. The performed analysis suggests that an improvement in the code should al− low adaptable calibration during the course of the eruption in order to mimic different emplacement styles in different phases. results in a driving force having an azimuthal compo− nent, and by following this component point by point a steepest descent path (SDP) is obtained [Favalli et al., 2012a]. This path can be derived through standard tools embedded in GIS software [e.g. Neteler et al., 2012], and an SDP−based analysis can forecast flow paths and drainage networks of subaerial water flows [Tarboton, 1997]. Lava flows are more complex than water flows because some of the originally flowing lava solidifies and modifies the pre−existing topography [Tarquini, 2017]. As an example, the DOWNFLOW code for the simulation of lava flows [Favalli et al., 2005] accounts for the feedback effect due to lava solidification by in− troducing a stochastic perturbation of the topography. The ways in which pre−existing topography controls the propagation of lava flows, and how this can be best implemented in lava flow codes, remains a subject of study. On the one hand it has been argued that the sub− tle details of topography are less important to include in a simulation if the lava flow emplacement is perturb− ing the topography [Wright et al., 2008]. On the other hand, there is a belief that technological and conceptual improvements in both lava flow modeling and topo− graphic surveys could reduce the simulation of the em− placement of lava flows to a deterministic process in which the tiniest topographical changes could be uni− vocally accounted for [e.g. Gregg, 2017]. In the present contribution, we report about a series of simulations started during the development of the probabilistic code for the simulation of lava flows named “MrLavaLoba” [de’ Michieli Vitturi and Tarquini, 2018], focusing on the effect of the pre−emplacement topog− raphy and on syn−emplacement topographic modifica− tions in the simulation output. We began our tests im− mediately before the beginning and during the initial phase of the 2014−2015 Holuhraun eruption and we completed our tests post−eruption as detailed field data about the actual lava deposit become available for com− parison [Pedersen et al., 2017]. Different pre−emplace− ment topographies with increasing quality were made available during the progression of the study. 2. THE BÁRÐARBUNGA VOLCANO AND THE 2014-2015 HOLUHRAUN ERUPTION The Bárðarbunga−Veiðivötn volcanic system consists of a glacier−covered central volcano and a 190−km− long fissure swarm [Einarsson and Sæmundsson, 1987; Johannesson and Sæmundsson, 1998; Thordarson and Larsen, 2007; Hjartardóttir et al., 2015; Larsen et al., 2015]. The central volcano is located within the north− western part of the Vatnajökull ice cap (Figure 1A), and the eruption frequency of the system in Holocene times is ~5 eruptions per century (Óladóttir et al., 2011). Re− cent lava flow fields formed in the Holuhraun floodplain in 1797 and 1867 [Hartley and Thordarson, 2013]. The 2014–2015 Holuhraun eruption was preceded by an in− tense seismic swarm caused by an impressive propaga− tion of a dyke from the central volcano to the northeast for ~48 km in about 2 weeks [Sigmundsson et al., 2014; Ágústdóttir et al., 2016]. On August 29, an initial 4−h long lava outflow was recorded along the old cone row, about 6 km to the north of the Dyngjujökull glacier. After a 24 h hiatus, the eruption resumed from the same eruptive fissure, and was continuously active until February 27, 2015 (IMO, 2015). During the volcanic unrest, a 65 m deep gradual caldera subsidence took place at the Bárðar− bunga central volcano [Gudmundsson et al., 2016]. The eruption was characterized by an average effusion rate of ~90 m3s−1, with an initial maximum of ~350 m3s−1. By the end of December, the effusion rate dropped to ~60–70 m3s−1, and continued to decrease from mid− January until the end of the eruption [Coppola et al., 2017]. The total bulk volume of erupted lava has been quantified in ~1.42 ± 0.06 km3 based on TanDEM−X DEMs analysis of pre and post−eruption topography (Dirscherl and Rossi 2018). There were prodigious emis− sions of gases during the eruption [Gislason et al., 2015; Pfeffer et al., 2018] and the cooling lava field contin− ued to outgas for months after the end of lava effusion [Simmons et al., 2017]. During the 6 months of continuous effusion, the emplacement dynamics changed significantly, and the flow field evolution has been divided into three main phases based on the dominant type of transport system delivering lava to the active flow fronts [Figure 1B, Ped− ersen et al., 2017]: phase 1, with channel−fed lava transport (August – mid−October); phase 2, with lava ponding (mid−October – November); and phase 3, with tube−fed lava transport (from December to the end of the eruption). The surface morphology of active lavas also changed repeatedly within a continuous spectrum from pahoehoe to ‘a’ā (in both directions, and also in a cyclic way). The analysis of the emplacement dynamic discussed in Pedersen et al. [2017] highlights that the lava supply was initially distributed through 8 large TARQUINI ET AL.


INTRODUCTION
Lava flows can pose a significant threat to human life and infrastructure, and forecasting areas which may be inundated by advancing lava is a critical issue [e.g. Kereszturi et al., 2014;Cappello et al., 2016]. Several lava flow simulation codes have been, and continue to be, developed [e.g. Cordonnier et al., 2016;Dietterich et al., 2017] and the scientific community is developing and promoting best practices to support further ad− vances and to better connect researchers and commu− nities vulnerable to lava flows [Harris, 2015;Harris et al., 2016;Mossoux et al., 2016].
During effusive eruptions, lava is erupted as a com− plex multiphase, multicomponent system made of a silicate melt plus a variable amount of solid fraction (crystals) and gas (exsolved volatiles). This mixture, as a whole, behaves like a non−Newtonian fluid, and moves downhill under the effect of gravity [Harris, 2013]. The interaction between gravity and topography results in a driving force having an azimuthal compo− nent, and by following this component point by point a steepest descent path (SDP) is obtained . This path can be derived through standard tools embedded in GIS software [e.g. Neteler et al., 2012], and an SDP−based analysis can forecast flow paths and drainage networks of subaerial water flows [Tarboton, 1997]. Lava flows are more complex than water flows because some of the originally flowing lava solidifies and modifies the pre−existing topography [Tarquini, 2017]. As an example, the DOWNFLOW code for the simulation of lava flows [Favalli et al., 2005] accounts for the feedback effect due to lava solidification by in− troducing a stochastic perturbation of the topography.
The ways in which pre−existing topography controls the propagation of lava flows, and how this can be best implemented in lava flow codes, remains a subject of study. On the one hand it has been argued that the sub− tle details of topography are less important to include in a simulation if the lava flow emplacement is perturb− ing the topography [Wright et al., 2008]. On the other hand, there is a belief that technological and conceptual improvements in both lava flow modeling and topo− graphic surveys could reduce the simulation of the em− placement of lava flows to a deterministic process in which the tiniest topographical changes could be uni− vocally accounted for [e.g. Gregg, 2017].
In the present contribution, we report about a series of simulations started during the development of the probabilistic code for the simulation of lava flows named "MrLavaLoba" [de' Michieli Vitturi and Tarquini, 2018], focusing on the effect of the pre−emplacement topog− raphy and on syn−emplacement topographic modifica− tions in the simulation output. We began our tests im− mediately before the beginning and during the initial phase of the 2014−2015 Holuhraun eruption and we completed our tests post−eruption as detailed field data about the actual lava deposit become available for com− parison [Pedersen et al., 2017]. Different pre−emplace− ment topographies with increasing quality were made available during the progression of the study.

THE BÁRÐARBUNGA VOLCANO AND THE 2014-2015 HOLUHRAUN ERUPTION
The Bárðarbunga−Veiðivötn volcanic system consists of a glacier−covered central volcano and a 190−km− long fissure swarm [Einarsson and Saemundsson, 1987;Johannesson and Saemundsson, 1998; Thordarson and Larsen, 2007;Hjartardóttir et al., 2015;Larsen et al., 2015]. The central volcano is located within the north− western part of the Vatnajökull ice cap ( Figure 1A), and the eruption frequency of the system in Holocene times is ~5 eruptions per century (Óladóttir et al., 2011). Re− cent lava flow fields formed in the Holuhraun floodplain in 1797 and 1867 [Hartley and Thordarson, 2013]. The 2014-2015 Holuhraun eruption was preceded by an in− tense seismic swarm caused by an impressive propaga− tion of a dyke from the central volcano to the northeast for ~48 km in about 2 weeks Ágústdóttir et al., 2016].
On August 29, an initial 4−h long lava outflow was recorded along the old cone row, about 6 km to the north of the Dyngjujökull glacier. After a 24 h hiatus, the eruption resumed from the same eruptive fissure, and was continuously active until February 27, 2015(IMO, 2015. During the volcanic unrest, a 65 m deep gradual caldera subsidence took place at the Bárðar− bunga central volcano [Gudmundsson et al., 2016]. The eruption was characterized by an average effusion rate of ~90 m 3 s −1 , with an initial maximum of ~350 m 3 s −1 . By the end of December, the effusion rate dropped to ~60-70 m 3 s −1 , and continued to decrease from mid− January until the end of the eruption . The total bulk volume of erupted lava has been quantified in ~1.42 ± 0.06 km 3 based on TanDEM−X DEMs analysis of pre and post−eruption topography (Dirscherl and Rossi 2018). There were prodigious emis− sions of gases during the eruption [Gislason et al., 2015;Pfeffer et al., 2018] and the cooling lava field contin− ued to outgas for months after the end of lava effusion [Simmons et al., 2017].
During the 6 months of continuous effusion, the emplacement dynamics changed significantly, and the flow field evolution has been divided into three main phases based on the dominant type of transport system delivering lava to the active flow fronts [Figure 1B, Ped− ersen et al., 2017]: phase 1, with channel−fed lava transport (August -mid−October); phase 2, with lava ponding (mid−October -November); and phase 3, with tube−fed lava transport (from December to the end of the eruption). The surface morphology of active lavas also changed repeatedly within a continuous spectrum from pahoehoe to 'a'ā (in both directions, and also in a cyclic way). The analysis of the emplacement dynamic discussed in Pedersen et al. [2017] highlights that the lava supply was initially distributed through 8 large TARQUINI ET AL. channels, which were active sequentially one at a time. As the emplacement progressed, resurfacing became important, as well as inflation and distal breakouts. Towards the end of the eruption, long tube systems, which seemingly resumed the pathways of the initial channels, played an important role.

THE MRLAVALOBA CODE
The Holuhraun eruption began in August 2014, when only an embryonic version of the MrLavaLoba code ex− isted [de' Michieli Vitturi and Tarquini 2018]. On the one hand, this was an evident drawback, because develop− ing the tool for performing the simulations during an ongoing crisis is not recommended [Harris et al., 2016]. On the other hand, the unusual characteristic of this eruption and the limited quality of the available topog− raphy at the onset of lava effusion (described later) were quite challenging, and stimulated the implemen− tation of specific modeling strategies.
In MrLavaLoba, erupted lava is itemized in parcels having an elliptical shape and (optionally) prescribed volume. Given the elliptical shape and the way a sim− ulation proceeds (explained below), we named our com− putational parcels "lobes", but these lobes are purely computational and probabilistic items, and must not be confused (nor compared) with real lava lobes. In the model, these lava parcels/lobes are intended to represent a volume of lava settled in a given position on the to− pography. Solid lava will always deposit along the prop− agation of the flow, which is controlled by the local slope. The code, however, does not directly account for any movement of lava parcels, which, in practice, are considered only when they have become solid. The simulation proceeds through the budding of new parcels from existing ones in a specific position determined by a probabilistic law influenced by the local steepest slope direction and by tunable input settings. In contrast with other probabilistic codes [e.g. DOWNFLOW, Favalli et al., 2005], MrLavaLoba can account for the volume of lava erupted. The code is intended to forecast the most prob− able area that may be inundated and the thickness dis− tribution of the lava deposit. One limitation of MrLavaL− oba is that, by neglecting the actual dynamics of the lava transport, the velocity of the progression downhill of the flow field is not considered, i.e. it is not possible with MrLavaLoba to forecast how long it may take un− til an area is inundated.
In the present work, beside setting a constant area and thickness for all parcels/lobes (hence also a constant volume), we have pre−defined the cumulative lava vol− ume in the simulations, and we account for thickness of the simulated deposit as a modification of the compu− tational domain. In this way, simulation results should better deal with syn−emplacement topographic modifi− cations. Crucial tuning parameters in MrLavaLoba, both varying between 0 − 1, are the parameters "lobe_exponent" and "α max " [de' Michieli Vitturi and Tarquini, 2018]. By setting lobe_exponent to 0, the budding is al− ways from the last/youngest lobe, and the progression of the simulation forms a single chain [no bifurcation is allowed, as in simple steepest descent paths in DOWN− FLOW − Favalli et al., 2005]. As lobe_exponent in− creases, instead, the budding becomes scattered throughout all the existing lobes, and a large amount of branching can occur, favoring flow spreading. The pa− rameter α max controls the extent of a random deviation from the steepest descent path. With α max = 1, the di− rection of a new lobe is the steepest descent path, and as this parameter decreases to 0, all directions have an equal probability, i.e. randomness is driving the direc− tion and the topography is not having an impact. There are additional tuning parameters the user can adjust, such as fixing the number of lobes per chain (discussed in more depth ahead). In this case, the cumulative num− ber of lobes are emplaced by iteratively creating many chains whose volume will sum to the prescribed cumu− lative volume. This setting can allow the user to con− trol the expected maximum runout, which is approxi− mated by the average length of the lobes times the number of lobes in each chain (all the other settings be− ing equal). In this case, the spreading can be controlled by the parameter α max . Thus, at fixed total volume and runout length, the spreading controls the area covered and impacts the thickness of the deposit (the greater the area, the thinner the average thickness).
At the onset of an eruption one cannot know the runout length, the deposit thicknesses and the cumula− tive volume. Studies of deposits from prior eruptions can provide values that have been produced from a given volcano or a suitable analogue, so it can be possible to explore a plausible spectrum of possible scenarios. This approach can be dynamically adapted during an erup− tion, to incorporate observations of the lava character− istics as they are made and communicated to the per− son running the MrLavaLoba model.

THE AVAILABLE DEMS OF THE HOLUHRAUN FLOOD-PLAIN
Four different DEMs were tested for lava flow sim− ulations at Holuhraun (Table 1). When the unrest began in August 2014, a 25 m−resolution DEM of the Holuhraun area was available to us. This DEM, here called "ÍSOR DEM", covers the whole country and is based largely on contour lines from 1:50,000 maps of Iceland but in some places it is based on more accurate data made by Icelandic Geosurvey (ÍSOR) [Víkingsson, 2008]. In the ÍSOR DEM, the floodplain north of Dyn− gjujökull was rendered as a completely flat area, in which the vertical gradient is captured in sparse, 2 m− height steps (Figure 2A). A local null gradient constitutes a significant problem when the purpose is to model a gravity−driven flow. To overcome this, the original contour lines were automatically re−derived along the steps, and processed using the DEST algorithm [Favalli and Pareschi, 2004;Tarquini et al., 2007]. The resulting improved DEM, here called "L−DEST DEM", rendered in the form of a 25 m−cell size grid, ( Figure 2B), largely re− duced the extent of completely flat areas (Figure 2A−B). A different 20 m−resolution DEM, mainly derived from remote sensing data (CARTOSAT satellite stereo images and COSMO−SkyMed radar images), usually called the "IPA DEM", referring to the IPA (Instrument for Pre−Ac− cession Assistance) programme of the European Union was available. The IPA DEM was tested ( Figure 1C), but the obtained simulation outputs (see section 3.1) sug− gested a lower accuracy compared with the ÍSOR DEM. Similarly, we tested the ASTER DEM [Hirano et al., 2003], obtaining unsatisfactory results.
A new TanDEM−X DEM from DLR (resampled at ~20 m resolution) called here "DLR DEM" ( Figure 1D) was available to us from November 2014 onward. The DLR DEM was obtained from scenes acquired on Febru− ary 2011 [Dirscherl and Rossi, 2018], and is assumed to adequately represent the pre−emplacement topography. However, the floodplain of Holuhraun is described in Pedersen et al. [2017] as covered by fluvial deposits from the Jökulsá á Fjöllum river, featuring decimeter to me− ter scale banks, bars and river terraces, which may have relief up to 7 m close to the central parts of the river channel. These fine−scale morphological features that can affect lava flow propagation would not be ade− quately preserved in a 20 m resolution DEM (discussed later on). We computed the elevation residuals of the other DEMs with respect to the DLR DEM, assuming it as ground truth ( Figure 3). The residuals show that the L−DEST DEM is more similar to the DLR DEM than the IPA DEM. The correction introduced by the DEST algo− rithm improved the similarity between the ÍSOR DEM and the DLR DEM by 11%, and this correction was enough to introduce a substantial improvement of the simulation results.

EARLY LAVA FLOW SIMULATIONS DURING DYKE PROPAGATION
During the 2 weeks−long north−eastward propaga− tion of the dyke, prior to the onset of the effusive erup− tion, the VORIS code [Felpeto et al., 2007] and an early version of the MrLavaLoba code were used to test pos− sible directions of lava flow advancement from three (arbitrarily chosen) potential vents located to the North of the Dyngjujökull glacier along the direction of prop− agation of the earthquake swarm ( Figure 4). Starting on August 31, in the initial hours of the eruption, low vis− cosity lava propagated almost radially from an old eruptive fissure as a series of fast, overlapping sheet flows. The first 24 h of lava spreading was reminiscent of our early simulations ( Figure 4), but the actual flow maintained this style of propagation only for a limited time, at a much smaller scale than the forecasted spreading in Figure 4.
By the 2 nd day of the eruption, the flow had already  self−organized in a channel which advanced steadily along a rather straight path, in spite of the very low to− pographic gradient (~0.2 degrees along most of the ini− tial 17 km of progression, Figure 5) and in contrast with the almost isotropic spreading forecasted by the preliminary simulations ( Figure 4). This behavior sug− gested that the local very low slope exerted a very strong control on the direction of propagation of lava channels [Favalli et al., 2009a]. Over the first few days, the fissure concentrated into three main vents concentrating the lava outflow to a few points. Thus, the initial simulations were carried out considering three equally supplied computational vents (labelled in Figure 5). Thanks to the effort of the "field team" at the University of Iceland [IES, Pedersen et al., 2017], observations and maps of the advancing flow and es− timations of the erupted volume were available to constrain the model inputs and to check the simula− tion outputs on a regular basis. Independent remote sensing effusion rate estimates were available from the MIROVA system [Coppola et al., 2016. As previously described, the ÍSOR DEM of the Holuhraun area, the only DEM available at the onset of the eruption, had a severe problem due to the com− pletely flat areas. This artifact promoted an almost isotropic spreading (i.e. a radial propagation) due to the absence of an elevation gradient, in stark contrast with the observed steady progression of the channel along a consistent path for several days ( Figure 5). The ÍSOR DEM was later refined obtaining the L−DEST DEM, however some flat areas were still present. To overcome this, we introduced in MrLavaLoba an "in− ertial factor" (ruled by the parameter inertial_expo-nent) which, when computing the direction of a bud− ding lobe, allows for variability of the direction of propagation of the parent lobe. In this way, in the ab− sence of a topographic gradient, the flow continues to propagate approximately along the previously estab− lished direction.
Below, values of some key parameters used for simulations of the initial 0.5 km 3 of lava discharged are reported and briefly discussed (see also Table 2). We considered lava emission from three vents (red pinpoints in Figure 5). Lobe area was set to 4 × 10 3 m 2 (i.e. 10 times a 20 m cell size); lobe_exponent was set to 0, so that simple lobe−chains are produced; the number NL of lobes per lobe−chain was = 1500, with NL yielding the runout (the higher NL, the higher the runout); number of lobe−chains NC = 600; inertial_exponent > 0 (to overcome the problem of flat areas). The settings above result in a large number of lobe chains propagating outwards from the vent. Bundles of lobe chains trending along similar (gravity−driven) paths compare well with the steady advancement of lava channels in the first phase of flow field formation [Ped− ersen et al., 2017].
An additional critical point which was explored in detail during the early simulations of the Holuhraun test case, was how to deal with the volume of the emplaced lobes (i.e. computational lava parcels) during the sim− ulation. As mentioned above, during the emplacement there is a continuous transition of mass from the flow− ing (liquid) system to the static (solid) deposit [e.g. Tar− quini, 2017]. Depending on the specific dynamic at work, solidification can repeatedly resurface levees [e.g. Sparks et al., 1976;James et al., 2007], or otherwise can slowly (but continuously) increase the thickness of a solid crust which insulates the molten interior of an ac− tive sheet flow (Hon et al., 1994). In both the above ex− amples, topographical modifications may become im− portant, such as in the case of the lava channel described in Tarquini and de' Michieli Vitturi, [2014], in which the lava deposit constitute levees thickened up to ~20 m in 5 days. In general, topographic modifications often re− sult in the classic "topographic inversion" [e.g. Giaco− mini et al., 2009;Hamilton et al., 2013]. When the sim− ulation settings imply the settling of many lobe chains, if late chains are derived on a topography which is mod− ified by the settling of previous chains, then the paths of new chains will tend to increasingly diverge from the initial steepest descent path, which is typically quite similar to the initial flow path.

TARQUINI ET AL.
6 FIGURE 3. Elevation residuals between the DLR DEM (resampled to 20 m), taken as reference, and the L−DEST DEM (blue) and the IPA DEM (red). Plotted residuals have been computed over the area inundated during the first week of the effusion (dashed line area in Figure 2D).
However, in the examples of lava flows described above (in both channelized and sheet flows), the actual structure of the flow tends to keep an established path, inhibiting flow spreading for a given time and promot− ing thickening. In channelized flows, thickening is pro− moted by the growth of levees which can lead to the formation of perched flows [Patrick et al., 2011;Tarquini et al., 2012a]. In mature pahoehoe lobes and in sheets 7 LAVA FLOW SIMULATIONS AT HOLUHRAUN, ICELAND   Figure 8B is available at http://demichie.github.io/MrLavaLoba/, where the source python code for MrLavaLoba is also available.
flows, instead [e.g. Katz and Cashman, 2003;Walker, 2009], thickening is due to an endogenous style of em− placement, lava tubes develop along the initial gravity− driven paths, and flow propagation becomes less sen− sitive to further thickening.
Indeed, MrLavaLoba does not directly simulate the formation of levees nor the roofing over a flowing lava. But we implemented specific settings to account for the "buffer effect" exerted by the structure of the flow. This is done by tuning the parameter "thicken-ing_parameter" (varying between 0 − 1). When the to− pography is periodically updated to account for the modifications due to the settled lobes [de' Michieli Vitturi and Tarquini, 2018] the thickening_parameteris accounted for. A high thickening_parameter implies a lower effect of the topographic change, promoting thickening over spreading. The opposite is true for a low thickening_parameter. In simulating the Holuhraun flow, a relatively high thickening param− eter was selected to account for the significant thick− ening of the actual lava deposit after an initial phase of lengthening and widening [Pedersen et al., 2017].

SUBSEQUENT LAVA FLOW SIMULATIONS
After the initial simulations carried out on different DEMs ( Figure 5), further simulations were carried out using the preferred DLR DEM and a more mature ver− sion of MrLavaLoba. The improved version of the code is computationally more efficient, as the topography can be updated at every single step (i.e. at every new lobe) without affecting the computational time, and the steep− est descent path is derived according to a more effec− tive algorithm. The improved code attains the same runout with a much smaller number of lobes per lobe chain (Table 2). Figure 6 illustrates the details of later simulations carried out by considering a fixed number of lobes per chain and, in turn, about half of the erupted volume ( Figure 6E) and the total volume ( Figure 6F). Simulation outputs are compared with the morpholog− ical characteristics of the pre−emplacement topography ( Figure 6D) and with the real thickness of the lava de− posit ( Figure 6B). From November 2014 onwards, lava ponding, inflation and distal breakouts become impor− tant [Pedersen et al., 2017;Kolzenburg et al., 2018], and to better account for the actual distribution of lava, es− pecially in the southern portion of the flow field, it was necessary to introduce additional computational vents ( Figure 6A) near the southernmost "distributary center" and breakout indicated in Pedersen et al., [2017] (their Figure 2B). The location of ephemeral vents can be viewed as a "data assimilation" process which integrates the observations collected in the field within the simu− lations. In this case the strong propagation southwards during phase 2 suggested that the breakouts originat− TARQUINI ET AL.  Table 2. ing c2.2, c2.3 and c2.4 were able to drain a significant amount of lava, and to account for this we set an ephemeral vent near one of the southernmost breakout mapped in Pedersen et al. [2017]. There are several as− pects of the progression of the real emplacement (schematically summarized in Figure 6A) that are wor− thy of being considered together with the other maps. To begin with, the initial channel (c1.1 in Figure 6A) does not follow the braiding drainage network up to the major confinement D1 ( Figure 6D), but appears to fol− low (and is confined by) the remains of an ephemeral drainage channel (not always noticeable at the DEM res− olution used here) along a more southerly path. The confluence between this path and the mapped drainage network ( Figure 6C Figure 6A), approximately 4 months after the onset of the eruption [Pedersen et al., 2017]. About the same time, the confinement D1 was eventually overridden by two lobes (pinpointed in Figure 6B). The post−eruption el− evation profile derived from a post−eruption DEM (Tan− DEM−X from scene on 28 April 2015, see acknowl− edgements) along the section e 1 −e 2 [ Figure 6B, bottom right inset, modified after Kolzenburg et al., 2018] shows that the contact between channels c3.1 and c1.1 formed a topographic minimum which marks the former boundary of the flow field (i.e. after phase 1). Interest− ingly, the performed simulations failed to reproduce the northward confinement of c1.1, as well as the dip at the contact between c1.1 and c3.1. After the early emplacement of channel c1.1, chan− nel c1.2 formed from a distributary center along c1.1 about 4 km away from the eruptive fissure, and devel− oped along a branch of the drainage network which crosses the center of the floodplain (indicated by the magenta arrow in Figure 6C). The zoom in Z1 shows the confluence between this branch (confined by D2, Fig−  ure 6D) and the main branch of the Jökulsá á Fjöllum river (confined by D3). The frames of Z1 show that both the map of the actual lava deposit thickness and the simulation in Figure 6E reproduce the two main linea− ments (d1 and d2) which drove c1.2 into the Jökulsá á Fjöllum river. c1.2 re−joined c1.1 shortly afterwards, forming a kipuka, which was later submerged but re− mained quite evident as a low thickness area in both the map of the actual lava thickness and the simulations.
With the exception of D1 (discussed above), major morphological confinements (D2, D3, D4 and D5 in Figure 6D) consistently promoted the thickening of the real deposit, and this effect is well reproduced in our simulations. What is not satisfactory in simulation out− puts in Figure 6 is that the lava progression southwards is much lower than the real one, and that the thickness of the simulated lava deposit is too high in the distal sector of the flow (compare thicknesses in Figure 6B− 6F). Simulations in Figure 6 are obtained by setting a constant number of lobes per chain, which implies a uniform distribution of lava lobes for any distance down−flow, or a uniform distribution of lava volume for any distance down−flow ( Figure 7A), because lava lobes have uniform area, thickness and volume. It is gener− ally not granted that all the lobe chains attain consis− tently a similar distance from the vent, because a straight lobe chains extension could be hampered by an unfavorable topography, with lobes remaining trapped in loops within a short distance from the vent. Similarly, in given circumstances, the topography can promote lobe chains scattering in diverging directions . In the present case, however, lobe chains tend to progress steadily along unimpeded paths and tend to converge towards a topographic bottleneck (the tip of the floodplain to the West). The bottleneck effect combined with the volume settling distribution pre− scribed in Figure 7A imply that an equal amount of lava volume must be settled across a progressively narrower area at increasing distances from the vent, originating an excessive thickening of the simulated deposit towards the flow front ( Figure 6F). To overcome this issue, we changed the MrLavaLoba parameter determining (i) the length of the lobe chains, and (ii) the consequent dis− tribution of the lava deposit according to the down−flow distance (Figure 7), as explained in the following.
MrLavaLoba allows to set the probability density function (PDF) of lobe chains length according to sim− ple laws (continuous or discontinuous), or to define it as a generic beta distribution, allowing for a full spectrum of possible PDFs of lobe chains length [de' Michieli Vitturi and Tarquini, 2018]. Figure 7 shows (left−side charts) that the average length of lobe chains decreases as the PDF changes from a constant length (the "max− imum" in Figure 7A), to a uniform probability between ~40% of the maximum value and the maximum value ( Figure 7B), to a uniform probability between zero and the maximum value ( Figure 7C); conjugated right−side charts show that the same final volume is distributed quite differently according to the down−flow distance.
The simulations in Figure 8 are derived from the PDF of lobe chains length illustrated in Figure 7, and all the TARQUINI ET AL. other input parameters set as in the simulation in Figure 6F ( Table 2). Figure 8 shows that the overall lava coverage and the thickness distribution are sub− stantially improved with respect to the simulations in Figure 6, suggesting that the transport of lava in the real flow, as a whole, decreased steadily (perhaps lin− early) at increasing distances down−flow ( Figure 7B). Beside the tips highlighted for the simulations in Fig− ure 6 (which also hold for the new simulations), it is interesting to highlight further points on the basis of the new simulations and other maps. The zoomed frame ZA1 shows that in an apparently flat and homogeneous sector of the floodplain (ZA1.1), 1−2 km downhill from the eruptive fissure, the final thickness of the lava deposit is rather convolute (dashed arrow in ZA1.2). Further comparisons with the maps in Dirscherl and Rossi [2018], that are not shown here for conciseness suggest that this feature, clearly not reproduced in our simulations, was formed during the last phase of the eruption. The zoomed frame ZA2 shows that the lava ponding caused by the evident morphological confinement (D4 in Figure 6D) is correctly reproduced in our simulations, but the main flow directions suggested by the drainage (ZA2.1) and by the final lava deposit (ZA2.2, dashed arrows) are not. In our simulations this morphologi− cal obstacle appears to originate a moderate "shadow effect" which limits the propagation of lobe chains according to given directions (ZA2.3). The zoomed frame ZB1 further highlights (i) the inflation of chan− nel c1.1 during phase 3, and (ii) the concurrent formation of a low lava deposit and low morphology track at the contact between c1.1 and c3.1 (red arrows in ZB1.3, see also elevation profile in Figure 6B bottom−right inset). The zoomed frame ZB2 shows the contact between channels c1.1 and c1.3, high− lighting that at the end of the initial phase a low lava deposit marked this contact (red arrows in ZB2.2, this feature is the same as the one in ZB1.3). However, ZB2.3 shows that, during the last phase, the lava transport−deposition was concentrated exactly along the same track, as if a syn−emplacement−formed low topography trap had captured the lava transport− deposition system. The zoomed frame ZB3 shows the same feature of ZB2, but shows data from Kolzenburg et al. [2018]. It is worth reminding the contrasting role played by former lava transport systems in ZB1 and ZB2/3 during the last phase of the eruption. In ZB1 a lava flow resumed endogenously along a lava transport−deposition system established several months before, while in ZB2/3 a new (portion of) lava transport−deposition system formed mainly by using former transport systems as bare topography. Although a better fit is obtained in Figure 8, per− haps the most difficult feature to be reproduced in our simulations was the lava progression towards the south during phase 2 (Figure 1 and 6A). Indeed, the real propagation of several channels during phase 2 (e.g c2.2−c2.4) is quite at odds with respect to the local direction of the drainage network. We illustrated above that MrLavaLoba can enhance the possibility of flow spreading against a low topographic gradient [de' Michieli Vitturi and Tarquini, 2018]. When we tuned the code to fit the southward spreading, we simultaneously obtained higher spreading towards all the other directions, resulting in a much poorer agreement between the real lava deposit and the sim− ulation output. It appears that different sectors of the flow field (in both space and time) were emplaced according to different dynamics, which can be repro− duced by the model only by utilizing different tunings of the code, while the code is currently writ− ten to allow for one set of parameters throughout one simulation.

DISCUSSION
The development of the code during the eruption at Holuhraun required particular care in addressing both (i) the effect of areas with very low slope (~0.2 degrees on average along the initial 15 km), and (ii) the impact of the syn−eruptive topographic modifications. It is inter− esting to further highlight the relation between the ac− tual dynamics at work [documented in the literature, e.g. Coppola et al 2017;Dirscherl and Rossi 2018;Pedersen et al., 2017 andAufaristama et al., 2018] and the mod− eling strategies implemented in MrLavaLoba to repro− duce lava emplacements influenced by these dynamics.
The existing literature suggests that, when the em− placement of a lava flow occurs on the relatively steep flanks of a volcano, the details of the topography are not crucial. At Mt Etna, for example, simulating lava flows on a DEM having a root mean square error (RMSE) in elevation of ~2 m [i.e. the TINITALY DEM,  is viable, and the output of a given simula− tion does not change significantly even if the accuracy of the DEM is further degraded [Favalli et al., 2009a; Tar− quini and Favalli, 2016]. The latter point is confirmed by Wright et al. [2008], which simulated the Mt Etna 1991− 1993 eruption based on a DEM having an accuracy sig− nificantly lower than that of the TINITALY DEM [~4 m, Wright et al., 2008]. In other volcanoes such as Mt Cameroon and Nyiragongo (DRC), the 90 m−cell size SRTM DEM [Rabus et al., 2003] was found to be good enough to perform reliable lava flow simulations Favalli et al., 2009b, respectively], de− spite its elevation accuracy is relatively low (e.g. it has an RMSE of ~8 m at Mt Etna, Tarquini et al., 2012b]. On the basis of our analysis of the available DEMs of the Holuhraun area (Figures 2,3), we infer that on the almost flat floodplain at Holuhraun, a simulation of the ongo− ing eruption using the IPA DEM (RMSE in elevation = 5.13 m, Figure 3) results in an essentially biased output (Figure 5a), while the same simulation (i.e. same input set− tings) using the L−DEST DEM (RMSE = 1.23 m, Figure 3) provides a similar results compared with the simulation using the reference DLR DEM (Figure 5b,c). However, im− portant differences are still present. The above suggests a much higher sensitivity of simulation output to DEM accuracy in the almost flat floodplain at Holuhraun. This outcome is not surprising, because where the elevation gradient is very low, a small perturbation of the DEM (i.e. a small error in el− evation) can have a very large effect on DEM−derived steepest descent paths. In practice, in flat areas a small error in elevation can determine the whole elevation gradient. On steep gradients, in contrast, the same per− turbation can be negligible, because the existing topo− graphic gradient dwarfs the contribution to the gradi− ent due to the perturbation in elevation.
The other interesting point is how to account for the observed emplacement style by tuning the input pa− rameters in the MrLavaLoba code (when possible). As highlighted above, one of the key points is dealing with the syn−emplacement topographic modifications. When considering entirely different effusive episodes spaced in time, the existing literature suggests that recent lava de− posits can act as a "repeller" for the path of new lava flows, as a result of the "topographic inversion" rule [Tarquini, 2017]. This is especially true when thick lava deposits form [Mt Etna examples are illustrated in Tar− quini and Favalli, 2010], and there is no specific topo− graphic confinement [Pinkerton and Wilson, 1994;Favalli et al., 2012a]. According to this, if the purpose would be to simulate discrete events, the settings within MrLavaL− oba should periodically account for the total change in topography. This can be attained by setting thicken-ing_parameter = 0, promoting spreading. But during a long−lasting eruption, things can work differently, and the Holuhraun eruption case is a good example of this. In the initial two phases (i.e. up to November 2014) the flow experienced an important lengthening followed by spreading with only minor resurfacing and thickening (which would suggest setting thickening_parameter to low values). In the last phase (from December 2014 to the end of the eruption), instead, the dominant processes were thickening and resurfacing in the middle of the flow field [Pedersen et al., 2017], and the latter emplacement condition can be attained by setting thickening_parameter to high values (i.e. thickening is promoted). In par− ticular, it seems that in the last three months of effusive activity, most of the newly erupted lava resumed to flow endogenously along the track of the early channels em− placed during September 2014. A similar behavior would appear impossible to recreate with the model if the lava transport system would have developed over a topogra− phy modified by the previously emplaced lava. To mimic a similar behavior with MrLavaLoba would imply setting thickening_parameter to 1. At the same time, however, a significant widening of the flow field occurred to− wards the North, and this part of the late progression of the flow field would be promoted by a relatively low value of the thickening_parameter, in contrast with the pre− viously addressed thickening in the middle of the flow field. In conclusion, the coexistence of both thickening and widening processes, points towards adaptable settings throughout a simulation (in both space and time, in case).
The analysis presented in Figures 6 and 8 further   13 LAVA FLOW SIMULATIONS AT HOLUHRAUN, ICELAND FIGURE 9. Extrapolation of the cumulative erupted volume based on the initial 5 weeks of observations. Data from the MIROVA system [Coppola et al., 2016. See also main text.
highlights the very close interplay between local slope (i.e. the existing drainage network) and lava flows. In general, during the Holuhraun eruption, flowing lavas tended to arrange themselves in channels trending along the ex− isting hydrological drainage network, as previously doc− umented in other cases of lava flow propagation in flat areas [e.g. Keszthelyi et al., 2004;Bernardi et al., 2015]. However, it is important to highlight the following ex− ceptions. First, the path and confinement (to the north) of channel c1.1 suggests that the direction of propagation of a lava flow can be the outcome of minor topographic confinements [Pedersen et al., 2017], unnoticeable at the resolution of the DEMs used here as computational do− main for simulations. In this case, the structure of the es− tablished lava channel, later on, can amplify the original minor confinement as a result of the self−confinement exerted by the levees [Tarquini, 2017]. Secondly, lava channels, in case of a very low topographic gradient, can locally propagate at odds with respect to the local slope (phase 2 Figure 6 A, C). The analysis of the volume settling distribution ac− cording to the down−flow distance (Figures 7,8) high− lights that this parameter can be crucial in order to bet− ter fit the real coverage and thickness of the lava deposit, suggesting a strong link between lava transport and de− position processes. In lava flow modeling, the direction of flow propagation has been successfully tackled through probabilistic approaches [e.g. Macedonio et al., 1990;Favalli et al., 2005;Felpeto et al., 2007], and this solu− tion has been often combined with other deterministic approaches covering the expected runout of flow units [e.g. Harris and Rowland, 2001;Chevrel et al., 2018] to obtain an improved, integrated tool [Wright et al., 2008;Harris et al., 2011;Mossoux et al., 2016]. Here we showed that, in addition to the direction and maximum runout, the distribution of the lava deposit (i.e. the p.d.f. of vol− ume settling) along the path plays a critical role ( Figures  6,8). We found that a volume distribution decreasing linearly from the vent up to the maximum runout pro− duces a much better fit than a constant volume distribu− tion throughout the path.
As a last point, it is worth highlighting that remote sensing data, in combination with modelling tools, can be used not only for near real−time hazard assessment [Wright et al., 2008;Vicari et al., 2011;Ganci et al., 2012;Cappello et al., 2016], but also to forecast the final emplacement during an ongoing event. The Holuhraun eruption showed that, if the driving factor has been elu− cidated (in Holuhraun, the collapsing caldera was pre− dictive of the lava effusion rate), the whole trend of the supply rate can be extrapolated on the basis of early ob− servations combined to a suitable model [Coppola et al., 2017, Figure 9], leading to beneficial forecasts of the fi− nal volume of the erupted lava. Remote sensing data and field observation could help in tuning the code in real time, according to different settings in different phases, in order to account for the different style of emplacement [i.e. channel−fed or tube−fed; Coppola et al., this volume].

CONCLUSION
We illustrated how the implementation of the MrLavaLoba code during the 2014−2015 Holuhraun lava flow emplacement promoted the development of specific modeling strategies. The general control exerted by the pre−emplacement topography on the direction of prop− agation of lava flows was already assessed by Favalli et al. [2009a] for the steep flanks of Mt Etna. Here we as− sessed the interplay between topography and lava flow in the case of emplacement on a very flat area, and we high− light the following points: 1) As a general rule, the lower the average slope of the pre−emplacement topography, the higher the impact of DEM accuracy in the simulation of lava flows. 2) The development of different types of lava trans− port systems (e.g. open channels or lava tubes) im− plies a different role of the syn−emplacement to− pographic modifications in the determination of the direction of propagation of the flow. Such different behavior can be mimicked (to a given ex− tent) in MrLavaLoba through optimized tunings. 3) The above point (2) further suggests that, to ac− count for changes in the emplacement style dur− ing an ongoing eruption, it should be necessary to change the input parameters during the course of an eruption as new observations are available. 4) Small scale morphological features (i.e. below the resolution of the DEM) can have a strong impact in the propagation of lava flow units, with impli− cation in the overall growth of the lava flow field. 5) The p.d.f. for the volume settling down−flow has a strong impact in the spreading of the flow and must be tuned in order to fit the real lava deposit. 6) The performed analysis suggests the introduction of further modifications of the code in order to ac− count for different styles of emplacement (i.e. dif− ferent tunings) during the course of the eruption.