Attempt to Model Lava Flow Faster Than Real Time: An Example of La Palma Using VolcFlow

: The eruption of Cumbre Vieja (also known as Tajogaite volcano, 19 September–13 December 2021, Spain) is an example of successful emergency management. The lessons learnt are


Introduction
Volcanic eruptions are one of the most disruptive processes on Earth. Their potential to transform the landscape and effects on public health make volcanic eruptions a complex process to manage. Volcanic eruptions are generally accompanied by a significant impact on the ecological and social environment, both local or global, and immediately or in the long term. Eruptions often force people living nearby to abandon their land due to the emplacement of lava flows [1][2][3][4][5] or by other volcanic products [6][7][8][9]. Volcanic explosions may fully obliterate an area (such as the well-known cases of the Krakatoa or Tonga volcanoes). Other than lava flows and explosions, eruptions may create different types of hazards [10,11]. Lahars and jökulhaups result from snow or ice melting [12][13][14][15]; toxic gas can come from lava degassing, lava interacting with a body of water, or from the fires started by lava [16]. Terrain collapses: earthquakes and landslides are also on the long list of processes an eruption may trigger [17][18][19][20]. Unlike other natural phenomena, volcanic eruptions may not be continuous; hence, they are frequently described in episodes with various processes and products within an indeterminate time frame of variable magnitude or intensity. Therefore, volcanic eruptions in populated areas require preparing, planning, and anticipating responses based on different hazard scenarios [21,22]. In the event of a volcanic eruption, emergency management must ensure that systems and services are in place to provide rapid and effective assistance. Faster than real time modelling is amongst the tools that may enable precise forecasts of some processes; hence, they contribute to better emergency management.
Live measuring the hazards of lava flows is a significant logistical challenge. Considerable time is required for processing, analysing, and interpreting the large amount of collected data. Nowadays synced-with-eruption lava flow assessment mainly relies on periodic surveys at different spatial scales [21]. Numerical modelling is usually carried out within pre/posteruptive scenarios either to provide hazard maps, for calibration of models, or inferring lava properties, as is explained below. The ecosystem of models is diverse, ranging from simplified to complex computer fluid dynamics to machine learning approaches. The outputs of such a wide variety of approaches are also reflected in the different types of products they provide and, consequently, the uses of such outputs. In the development of the temporal phases in which a high-impact lava flow event can be found [23], lava flow modelling plays a relevant role: (i) Before the crisis models are calibrated and validated based on past events [24][25][26][27]: Their outputs are useful to inform and educate about hazards and risks [28][29][30][31][32][33][34], and they provide knowledge to support measures to face the most likely scenarios, from expected volcanic activity to its spatiotemporal probability of emission points [35], as has been conducted since the 1980s [36]. (ii) During the eruption: Models can be used to predict the path and emplacement of lava flow [25,[37][38][39][40] and to assess potential scenarios or mitigation measures (such as to deviate the lava flow [41][42][43]). (iii) After the eruption, hazard maps should be updated. Simulations require new topographic data due to the different kinds of landscape changes. Several studies have been carried out in the Canary Archipelago to assess volcanic hazards [44][45][46][47][48], including on the Island of La Palma [49,50]. All published models are in the preparation phase, and few have been published in scientific journals to date [51][52][53]. Some have been lately released on the websites of IGN and La Laguna University (https://www.ign.es; https://volcan.lapalma.es/pages/visor, accessed on 22 November 2022). Their results show significant differences from the actual lava flow extent, showing a compromise between timely forecasts and predictability [21,54], but is also linked to the kind of outputs the models provide.
Our study was carried out to explore if the numerical simulation of lava flow can be used for forecasting during an emergency. The aim was to construct such an approach based on minimum preliminary tools and information so it may aid in identifying where are efforts most needed and what can be implemented. Because data specific to the recent La Palma volcanic event are not fully open-access, we relied solely on publicly available information.
We used a simplified model that is easily accessible and can be quickly analysed. The author of the model (Karim Kelfoun) kindly provided the original code files for this experiment. In this context, the study was performed as follows: (1) Evaluate the different approximations that could be obtained for the rheology of the lava from remotely obtained information. Such information included the evolution of the lava's extent during the first seven days of the event (19-26 September 2021). This time limit was set to avoid the modelling of the lava flowing down the coastal cliff. (2) Model the process of lava emplacement for the first five days, and add two more days as validation. We performed calibration using the best geometric fit assuming different scenarios and hypotheses. (3) Evaluate the different execution times of the proposed scenarios to determine the possibility of achieving a total calibration and simulation time shorter than that of the lava flow emplacement.

La Palma Island and Its 2021 Eruption
There is abundant literature about the geology of the Canary Archipelago, including that of La Palma island. Geological maps on a 1/25,000 scale and their reports are available from the CN-IGME website (http://www.igme.es, accessed on 22 November 2022; Figure 1). Although many issues regarding the geological origin of the archipelago are still far from being fully understood, their volcanic nature is beyond doubt. La Palma island comprises two clear separated regions. The northern part is mostly Pleistocene, whereas the southern part is covered by Holocene eruptions. It is generally assumed that this setting is also valid for the overall volcanic structure, which is mostly located underwater, raising from the seafloor at around −4000 m and reaching up to 2426 m above water at its highest peak. The morphology of the island is deeply marked by mega landslies [55] and an eruptive configuration compatible with that of ocean ridges [51]. Some authors remarked that there has been no eruptions in the northern part of the island for the last 125 ka, although they also recognize the difficulties of reconstruct ingthe island's eruptive history beyond 100 ka [51]. They have also proposed different settings to explain the geological evolution of the island. Historical eruptions are characteristic of ocean island ridges, involving lava flows and moderate explosivity, occasionally enhanced by the presence of evolved magma or due to the interaction of lava with water bodies [51]. In the historical record, at least eight eruptions (nine counting the one in 2021) have been identified, having an average duration of around 60 days [56] (63 considering the eruption in 2021). The most recent eruption started on 19 September 2021, at 15:12 h WEST. The 2021 La Palma eruption is classified as a fissure-type basaltic eruption, dominated by Strombolian activity with episodic phreatomagmatic pulses, similar to the preceding historical events [57]. The eruption opened several vents as it progressed, lasting up to 85 days. The numerous monitoring stations deployed in the Aridane Valley area by different national and local organisations (e.g., IGN, CN-IGME CSIC, and INVOLCAN) allowed the early detection of the eruption [11,51,58,59]. Monitoring such pre-eruptive signs was key to undertaking early measures to evacuating the many inhabitants of the area and bringing them to safety. Monitoring the volcanic eruption during its evolution was also essential to further providing first responders and managers the means to tackle the overall emergency, such as modifying safety perimeters; providing alerts for lava flows, gases, and ash; and organising other activities, such as the relocation of inhabitants or granting access to water supply for inhabitants, crops, and livestock. One such monitoring task included the daily update of the lava delineation paths and expert-based analysis of the most likely pathway to be followed by the lava [60]. Responders, scientists, and experts working during the catastrophe are still studying the event to better understand how the overall actions contributed to the safety and well-being of the inhabitants as well as means to improve the management of future crises.
This event in 2021 differs from previous historical ones in the number of people that were directly affected (more than 2300 according to different press releases). From an economical perspective, this eruption was the largest in terms of damage in the Spanish history, as well as the first time the public insurance sector compensated for volcanic damage (personal communication from the Managing Board of the Consorcio de Compensación de Seguros). New eruptions are expected, considering the island is located near an active volcanic. Hence, any tool to further support first responders might also be useful for such future scenarios.

Materials: Geospatial Information
In the course of an emergency, the information available to determine the evolution of the lava emplacement can be classified into two types: pre-event and syn-eruption. The pre-event spatial information includes the topography and its derived morphometry (slopes, curvature, and slope orientation) in the form of a numerical matrix (or digital terrain models (DTMs)). This terrain information is eventually altered by the emplacement of successive lava flows along the course of the effusive process. It is also necessary to take into account human infrastructures (houses, roads, and bridges) in which the heights and strength of the buildings or infrastructure determine their role as retaining walls against the progress of the lava. Syn-eruption information is collected through the use of airborne remote sensing systems or satellite platforms to locate the emission points of the volcanic products; to determine their emission rate; and to delimit the extent, velocity, flow, heights, temperatures, and rheological properties of the lava flow [61][62][63]. On-site sampling was also carried out by different scientific and emergency authorities, but data have not been fully released. The quality of the remote surveys is determined by the ground sampling distance (GSD) at which the image was taken, the off-nadir angle, and by the sensor type and the sensor's own precision and accuracy. Synthetic aperture radar (SAR), thermal, infrared (IR), and optical sensors are the most commonly used.
The IGN is a public body responsible for monitoring volcanoes across Spain, as well as for distributing basic geographic information and other duties. The IGN provides a free download web service (CNIG) for different types of data, including the DTMs of the natural ground (DEMs), buildings, and vegetation (DSMs). Within the available products, we used the newest digital terrain models (DTMs) derived from LiDAR airborne campaigns. Bare ground terrain can be downloaded with a 2 m pixel resolution (known as MDT02, hereafter DEM), and a 2.5 m pixel resolution DTM is available for buildings (known as MDSnE2.5). The orthometric height data of the IGN MDT02 raster have an altimetric accuracy 95% better than that of 0.2 m (RMSE Z, [64]). These models were obtained from the average elevation of the terrain class points sampled by the LiDAR flight for the National Aerial Orthophotography Plan (PNOA) with a density of 0.5 points per m 2 and with double quality control (in the LiDAR point cloud and in the derived raster). The IGN MDSnE2.5 raster represents the height relative to the ground of the building (LiDAR class 6). The full bare ground and buildings can be obtained by adding both (hereafter, DSM). The geodetic reference system for these data in the Canary Archipelago is REGCAN95 UTM 28N.
Geospatial syn-eruption image surveying is strongly affected by the prevailing meteorological conditions in the area and the density and development of the eruptive ash cloud. Therefore, satellite remote sensing systems were the first resource used to assess the extent of the eruption and its overall evolution. Additionally, after the fourth day of the eruption, the emergency management authority (PEVOLCA) in collaboration with the Cabildo Insular de La Palma, the CN-IGME and other public entities started the collection of detailed information using low-altitude UAVs to avoid exposure to the cloud cover, ash, and gases from the volcano.
The emergency management service (EMS) initially decided to use very-high-resolution (VHR) optical satellite images (Airbus DS SPOT-7, GSD 1.5 m). However, due to poor coverage, SAR-based imagery was used instead as a major source of remote data. One day after the onset of the eruption (20 September 2021), the Copernicus EMS rapid-mapping team (activation EMSR546, [60]) provided the interpretation of high-resolution SAR images. Such data were gathered by the COSMO-SkyMed Second Generation constellation (GSD 1 m/pixel, ©ASI, 2021) and by Copernicus Sentinel-1 A/B (GSD 10 m/pixel, EU ESA). The interpretation provided the activation extent Mmap (AEM): a delineation of the lava flow and the buildings, roads, and infrastructure that had been affected. As the event evolved, the AEM was updated at least daily. The interferometric processing of SAR images allowed the assessment of terrain deformations at any time of day or night including cloudy areas. Many other sensors and platforms were used when the cloud or dust coverage was very low, which was rare. These observations included infrared, optical, and multispectral from ASTER, Pléiades-1 A/B Landsat-8, Sentinel-2, and GeoEye-1. For the seven-day period analysed in this study, we used a set of 10 lava flow outlines from EMSR546 (Table 1, Figure 2a), which are still available from the [60] collection. The Cabildo Insular de La Palma contributed to the emergency response team implementing a mapping plan supported by UAVs, using optical mosaic images with centimetric resolution. Their products have been available since 24 September, but only four lava flow outlines were available in the catalogue (Figure 2b) for the time interval studied in this work up to 27 September. Out of these 14 lava outlines (Table 1), we measured the area [WL] i and the linear distance L Ai travelled by the lava from the vent ( Figure 3) at each moment t i . This set of 14 perimeters was jointly used to morphometrically assess the evolution and rheological properties, and to undertake the numerical calibration of the lava flow, as explained in the following sections.
The time frame of this study was limited from the beginning of the eruption (19 September 2021) until 27 September 2021 for two main reasons: the first days represent a pristine scenario, where the fewest assumptions are needed to start the numerical model; by 28 September, the lava reached the coastal cliff, forming a lava delta that requires modelling a significant boundary condition, which was outside of the scope of this study.
For the temporal sequencing of the lava emissions, the first days were considered as a point-source type on an elongated cone. This assumption seemed appropriate as the fissural eruptive vents were located very close to each other, in the Cabeza de Vaca area. The fissures followed the northwest-to-southeast direction, forming a string of small cones around the vents (Figure 2, V1 and V2 in Table 2). On 25 September, a new flow was produced after the collapse of the southwestern flank of the cone. A less viscous and more abundant lava followed through a corridor, running downhill in a westward direction ( [51,65] Figure 2 and V3 in Table 2), and ceasing the contribution of the initial emission points or rendering their output as negligible. This later source was also modelled as a point source (V3).
Calibrating and validating the simulation of the lava evolution can be improved using the distribution of the mass over the extent of the flow. This distribution can be quantified in the form of thickness or height of the lava layer. Unfortunately, the only available information on lava thickness was provided by Cabildo Insular de La Palma (https://www.opendatalapalma.es/, accessed on 21 February 2022. seven days after the start of the eruption, 25 September. Moreover, this information, although obtained by means of UAV with a very high resolution, was only partially open as an RGB-reclassified image of 2/3 of the lava flow.

Methods
A key objective of volcanology during an eruptive crisis is to provide forecasts of the emitted lava streams, including the time they take to reach populated areas, infrastructure, or any other area of interest. This information is extremely useful for emergency services in managing the evacuation of people and in mitigating losses [66]. In the last decades, several calculation and numerical simulation tools have been developed based on the thermorheological behaviour of lava in its evolution to reproduce its trajectory over the terrain's topography [21,[67][68][69]. The main challenge in designing these models is to simulate the complex thermo-fluid-mechanical interactions that determine the morphology, distribution, texture, thickness, and extent of the emplacement process of a lava flow.
Several of the most widely used computational models have been tested and validated in benchmarking exercises [67,68]. The numerical approach of mathematical models has become one of the most important tools for assessing volcanic hazards such as lahars, pyroclastic surges, ash clouds, and lava flows, amongst others. Specifically, the simulation is handled from different strategic perspectives for the evaluation of the evolution of lava flow emplacement. Some of the first computer codes were based on cellular automata schemes to solve the Navier-Stokes equations with a non-Newtonian rheology [38]. They later evolved into more complex algorithms such as MAGFLOW [39,70,71], SCIARA [72], FLOWFRONT [73,74], and MOLASSES [75]. Other models focused on a more complete thermorheological model over a one-dimensional channel, such as FLOWGO [76,77]. The knowledge gained in both lava physics and its mathematical and numerical modelling has enabled 3D simulations (LavaSIM [78]); the introduction of randomness as a way of handling uncertainty (DOWNFLOW [79]); or the description of the conservation of depth-averaged velocity, temperature, and flow thickness (VolcFlow [80,81], RHE-OLEF [54], and IMEX-SfloW2D [82]). Today's computer capabilities are allowing further advances, pushing the complexity of models (3D CFD toolbox of Open Field Operation and Manipulation (OpenFOAM) [83]) and their integration into commercial code (FLOW-3D [67,68]). Moreover, results are been presented of the integration of hydrodynamic models in particle-mesh-free methods (GPUSPH [84,85]; NB3D code [67]). With recent GPU-based modelling, computing speed has largely improved compared with that of CPU-based calculations [86].
Despite the long duration of the eruption of La Palma, compared with more recent historical eruptions on the island [57], the data generally available so far are either very limited, scarce, or are still in the process of being refined and studied. Therefore, a key objective in this study was to build a model in a comprehensive and expeditious manner based solely on readily available public-domain information. Our hypothesis was that the model would allow near or faster than real-time forecasting of the lava flow hazard during the days between 20 and 27 September 2021 of the Cabeza de Vaca eruption, within the constraints of the available information. To achieve this objective, we used a computationally fast approximation based on depth averaging, in which the molten fluid is considered as a continuous incompressible material of constant density moving over a complex topography. The following assumptions and simplifications were also required [87]: (i) vertical accelerations are negligible compared with the dynamics in the horizontal direction, (ii) the force balance in the vertical direction is expressed by a hydrostatic pressure balance, (iii) the rheology is simple enough to be estimated with the available data, so (iv) homogeneous flow and rheology properties are assumed throughout the vertical section and (v) isothermal 2D behaviour, and (vi) lava behaves as a single-phase molten material.

Inferring Lava Flow Rheology
The numerical model used must reproduce the dynamic behaviour of a viscous mass. This is governed by the relationships between the forces that cause the lava flow to move and those that slow it down due to cooling. Therefore, it is necessary to incorporate the amount of lava flowing into the system and its deformation behaviour under stress.
Walker [88] pointed out that the main factors controlling the extent and type of lava flow are its viscosity and effusivity rate, i.e., the rheological characteristics of the lava and the flow rate at which it is emitted. Numerous studies on lava [89] have shown the clear dependence of its characteristics on temperature, heat flux during cooling, dissolved volatiles, crystal concentration, and composition, although density shows little variation between lava types, typically ranging from 2100 to 2800 kg/m 3 [90][91][92] for basaltic lavas under atmospheric conditions, even though this value remains almost constant when the lava has solidified. The thermal conditions during cooling control both the change in rheological properties while the lava is flowing and the crystallinity of the solidified material after cooling.
How molten siliceous materials, which form lava, deform or flow in response to applied forces or stresses is not as simple as for Newtonian fluids with constant viscosity and is independent of any external stresses applied to it. Lava flow generally behaves like a rigid body when subjected to low stresses, but flows like a viscous fluid when subjected to higher stresses, i.e., like a Bingham fluid [93,94]. In field work on active lava [95,96], the behaviour of lava has been characterised as non-Newtonian flow. Measurements from these early works indicated that lava flow behaviour is pseudoplastic and can be approximated by a Bingham fluid. Phenocrystals and gas bubbles within the lava produce a stress threshold R 0 that must be exceeded for it to move. Once this elastic limit is exceeded (R 0 > 0), the flow can be characterised by an apparent viscosity. Conversely, a flow can be thin enough that the driving stress drops below R 0 , stopping the flow. The value of R 0 also influences the final thickness of the stopped lava. Although lava usually maintains its Bingham viscosity as the lava moves, retaining its characteristic aa or pahoehoe class, the yield stress and Bingham viscosities sometimes increase as the lava crystallizes and the flow velocity decreases, leading to the transition from pahoehoe to aa class, which seemed to have been observed in the first lava flows produced in the eruption of La Palma [51]. This can be quasi-sudden if the flow encounters an abruptly steepening change in slope, which causes a sudden increase in shear stress (steepening of the slope), thus driving the transition from pahoehoe to aa. However, this was not observed in La Palma during the study period, because the slopes of the terrain occupied by the different lava flows, from their point of emission, were almost constant or decreasing.
The viscosity in lava depends on many factors. The controlling factors are the silica content; the volume fraction of crystals; crystal size, shape, and, mainly, the temperature [97,98]. Basaltic lavas, at about 1000-1200 • C, such as those emitted by the volcano on La Palma (PEVOLCA Scientific Committee Report 20 September 2021), during the study period. The temperature of the lava remained almost constant during the first days (as shown in the thermal images of the lava taken with UAVs by the Cabildo Insular de La Palma and the IGME-CSIC (https://info.igme.es/eventos/Erupcion-volcanica-la-palma/videos, accessed on 22 November 2022)). The lava had low contents of crystals, silica, and dissolved gases, which gave them a relatively low viscosity compared with lavas of the dacitic or rhyolitic type [92]. Its apparent viscosity could be approximated by taking the crystalline fraction and its maximum value below which mass movement is allowed, using the Einstein-Roscoe relationship [99,100], which requires knowledge of this fraction and the temperature at which the lava flow is found, as well as other constants involved in the relationship. If the crystalline fraction is unknown, it can be estimated from the temperature and composition of the lava by taking the initial crystalline fraction (at the time of vent outflow), the total additional amount of crystallization occurring during the flow, and the solidification temperature [101].
However, although there are techniques that allow lava viscosity to be estimated in the field [102,103] or in the laboratory [104][105][106][107][108][109], these approaches require experimental analyses that, although accurate, are costly and time-consuming and would hinder reaching the real-time target. An alternative that is faster and commonly used in both local and planetary preliminary studies to estimate viscosity and yield strength is based on the morphological dimensions of lava flows. Initially, Nichols [110] proposed using Jeffreys' empirical formula [111]: which relates viscosity η 1 to lava height h and velocity u for a magma free of crystal content and providing a modified second formula: which incorporates the rate of lava emission Q from the source and a coefficient n that varies according to the width of the lava flow (n = 3 for wide flows or n = 4 for narrow flows). A simplified version of the former formulas [112,113] in which the influences of lateral flow spread and slope effects are omitted is: When lava is placed through a volume-limited flow [114], the effusion rate Q [115] is related to the length L, width W, and thickness h of the lava flow; the volume of its entire mass; and, for a given instant, is divided by the time elapsed since the start of the event. Therefore, the emitted volume Q per time unit is the variation in the thickness distribution h(x, y, t) at each point (x, y) within the spatial 3D region Ω inundated by the lava flow over time, can be approached by the perimeter-limited-area (WL) of the lava flow and divided by the elapsed time ∆ t since the eruption started and the time t at which the data were taken: However, when stopped by cooling-limited flow, thermal diffusivity is involved through the Grätz number G z , a dimensionless quantity used in volcanology [116] that expresses the balance between the heat transported by lava and the heat lost by thermal conduction. The value of G z is bounded between 100 and 300 when the lava flow is controlled by cooling. However, many flows are volume-limited and not cooling-limited, so the estimated values of G z in these cases are higher than expected. With these conditions, using expressions (5) and (6), assuming u is the average lava flow velocity, we propose making a preliminary estimate of the morphometry-based emission rates using: and: In practice, Jeffreys' expression, (1) or (2), is applied under laminar flow conditions, and when the lava flow sections have a width much larger than the depth or thickness of the flow [117] and assuming Newtonian behaviour of the fluid. Using one of the two, η 1 or η 2 gives an a priori estimate of the apparent viscosity [101] from depth-averaged measurements of the heterogeneous lava flow by taking an average thickness and velocity, an estimated value of its density, and the average slope of the terrain if it does not undergo major changes in relief. These values can be quickly obtained in the field and used in Equations (1)-(4).
The flow velocity v can be globally estimated as an average velocity from the interpretations of satellite images provided by emergency mapping services. Furthermore, the average velocity is obtained because, for small strain rates or low velocities, the yield stress can give rise to an apparent viscosity that is much higher than the actual viscosity. Thus, measurements of the geometry of the lava flow evolution, as in its expansion or in the main direction of advance (Table 1), allow an a priori estimation of the average velocities, of expansion u W Li = [WL] i /∆t 0i or average area expansion rate, and of the advance u i = L Ai /∆t 0i , with ∆t 0i = t i − t 0 being the time elapsed from the beginning of the eruption and the ith time (Table 1); and, likewise, of the instantaneous velocities of expansion ([WL] i+1 − [WL] i )/∆t i+1i , or instantaneous area expansion rate, and of advance (L Ai+1 − L Ai )/∆t i+1i . The volume values thus obtained, using the frequency distribution of lava thicknesses, can also be used as an approximation of the products h[WL] in Equations (4) and (6) to calculate the lava flow rate Q.
The precise area extent of the lava flow (WL) that is assessed at each time of the event using aerial or satellite interpretation to determine the lava outline, and the distribution of lava heights that require photogrammetric flight and DEM prior to the eruption, may be time-consuming or not feasible to obtain due to local visibility or safety conditions during real-time estimation. Hence, their values can be approximated with estimates based on the length L, average width W, and height h of the front of the lava flow. The value of the emitted volume, using the product of these quantities V i h[W L] i , at the ith time (Table 1), can be corrected with a shape factor [118] that reflects the difference in volume due to the geometrical irregularity of the lava flow. Considering the possible bias in the frequency distribution or histogram of lava thicknesses, the height h to be used in the product to determine V i is the mean <h> or the median h 50 , although, given the possible difference in results, it is possible to directly use the histogram of heights. If h * j is a characteristic point of each jth class, i.e., its midpoint, with f j frequency in the histogram, then the volume i , assuming that the lava thickness distribution is preserved over time. The viscosity value thus obtained disregards important physical and chemical processes that can affect rheology [119,120]. Lava flows show complex behaviour as viscosity varies with temperature, particularly due to the crystallisation induced by heat losses during lava emplacement. Viscosity will increase from its outlet to where it stops on several orders of magnitude. Therefore, lava viscosity is transient, meaning the validity of interpretations made from morphological methods must be considered with caution [121].
Once the above considerations have been given regarding the morphometric estimates of the advance velocity of the lava flow u and its flow rate Q for each time when information is available, it is now possible to estimate the viscosity over the observation time considered in this study. Using Equation (1), the thickness h is replaced in the numerator by its mean value <h> or the median h 50 , and the average velocity value u i , estimated at each i-th time, is applied in the denominator. The viscosity of Equations (2) and (3) can also be estimated taking the thickness h as its mean value <h> or the median h 50 , but using one of the flow rate estimates Q (Equations (4) to (6)) obtained in the denominator. However, considering the formulation used in this section for the previous estimates of velocity u i , volume V i , and flow rate Q, there is a functional relationship between them. So, for the estimation of the time evolution of the viscosity, only the rate Q obtained with Equation (5) was used because, if Equations (4) and (6) are used, the result is simplified, losing its variability over time. Despite these considerations, the morphological viscosity estimates for the La Palma eruption are a proxy that may prove useful for the purpose of real-time calculation and give a conservative estimate of its displacement and extent to inform an early warning management system. During the first days of the eruption in La Palma, the lava flow had little heat loss with an almost constant temperature throughout; the slope was almost the same along its westward path, except in the initial section near the point of emission; the density was constant; and the thickness, which remains relatively constant, can give rise to small variations in apparent viscosity, more due to structures that locally impede the passage in the path of the lava flow in favour of the topographic gradient than to internal rheological changes in the lava.

Numerical Approach and Simulation of Lava Flow
Simulations in this study were carried out with VolcFlow, a numerical code that was developed at the Laboratoire Magmas et Volcans, a joint research unit of Clermont-Auvergne University, affiliated with the Centre National de la Recherche Scientifique from France (CNRS, UMR 6524). The mathematical model is based on the depth-average approximation of shallow-water depth-averaged equations of mass and momentum balance, as is the case for other geophysical flows. VolcFlow is able to simulate various rheologies (plastic, viscous, Coulomb, etc.), and it is able to calculate the advection of additional variables (density and temperature) in modified versions. From its beginning to the present, Vol-cFlow has been successfully used for simulating debris avalanches [80,[122][123][124][125], pyroclastic flows [81,126], tsunamis [127], and lava flows [33,52,128]. Additionally, it has provided very accurate and fast results in benchmark tests [21,67,68], and it is very well documented that it is very well balanced in terms of complexity and usability. Therefore, it was considered as a good candidate, yet not the only one, to explore its use during a volcanic crisis for forecast purposes.
VolcFlow is a computational program that can be classified as deterministic simulation software developed to numerically simulate the emplacement of lava flow by modelling the physical principles of heat transfer and fluid mechanics [129]. The mathematical model incorporated in VolcFlow is based on the general depth-averaged mass and momentum equations or shallow-water equations (SWE), as hypothesised by Saint-Venant [87], assuming that the horizontal scale is much larger than the vertical scale. Using a topographiclinked coordinate system, with x and y parallel to the local ground surface and h vertical, the vertical component is neglected or discarded: where the subscripts in (8) and (9) denote the components in the X and Y directions [33,80], respectively. The SWE (7)-(9) are used to solve h = h(x, y, t), the thickness, and u = (u(x, y, t), v(x, y, t)) (the lava flow velocity) at a point of coordinates (x, y) at instant t, for a fluid with constant density that moves over a terrain of slope α(x, y) under some initial and boundary conditions. This system of partial differential equations forming the SWE are numerically solved by the finite difference with a Eulerian upwind integration scheme in MATLAB ® , a computing platform. This numerical scheme can handle shocks, rarefaction waves, and granular jumps and is stable even in complex topography and on numerically "wet" and "dry" surfaces. The isothermal version of VolcFlow used in this study allows simulating lava flows under conditions of constant Bingham rheology over time, independent of the cooling processes of the lava that could modify them. This rheological behaviour is incorporated in R = (R x ,R y ), which expresses the basal shear stress that depends on the rheological behavior model that is included in the second member of the SWE to reproduce the lava flow. In the VolcFlow program, the effect of several rheology models is implemented in the depth-averaged form as: (i) The Coulomb friction that relates the shear stress R to the normal stress at the base of the lava and to the basal friction angle between the lava and the original soil, which is formulated through the earth pressure coefficient k act/pass [130] or the ratio of the stress parallel to the ground to the normal stress of the ground: This lateral stress ratio, according to Mohr-Couloulomb theory, relates the depthaveraged horizontal normal stresses (diagonal components of the stress tensor) to the normal stress averaged in the vertical direction. The use of a scalar coefficient k act/pass pre-serves the invariance in the OXY plane, as well as the symmetry of the stresses. According to the sign of the numerator, in this coefficient, "−" is applied as active for divergent lava flows (k act ) and "+" is applied as passive for convergent flows (k pass ). When ϕ bed < ϕ int , the lateral normal stresses where the flow converges usually exceed those where the flow diverges by a factor of 2 to 10. The exception to this behaviour occurs if the bed has a maximum roughness; in this case, ϕ bed = ϕ int [130] and (10) simplifies to a single-valued expression for converging or diverging flow. In this case, the uniqueness of this value indicates that the material moves downslope as a slab, without thinning or thickening. (ii) The viscous friction is related to the dynamic viscosity η, with the speed being inversely proportional to the thickness. This is evident when a thicker casting moves slower than a thinner one under the same conditions. (iii) When the velocity is high, the Voellmy rheological model is used, in which the shear stress is proportional to the square of the velocity, through a coefficient ξ of Voellmy that represents the effect of turbulent dispersive pressure and/or collisions. (iv) A fluid is able to flow (i.e., deform indefinitely) only if it is submitted to a shear stress above some critical value or yield strength R 0 . Thus, a flowing Bingham lava is described when plastic and viscous behaviour act simultaneously. Bringing these models together in the shear stress R term in the SWE gives: With depth averaging and the resulting two-dimensional reduction in the originally three-dimensional problem, the CPU consumption times are low. However, height averaging implies assuming a constant rheology behaviour in the vertical profile, which limits the real characterisation.
To carry out a simulation with VolcFlow, the user supplies the program with data on the characteristics of the lava: (i) the rheological properties: internal friction angles ϕ int and basal ϕ bed , viscosity, and cohesion; (ii) density; and (iii) the volume emitted throughout the simulation time span. The DTM (DEM or DSM) on which the flow simulation is performed is specified as an ASCII file in the calculation script. During the simulation, at each time step, the computer code integrates the SWE system to obtain the distribution of thicknesses, velocities, and derived variables (pressures, stresses, mass, etc.) on the spatial discretisation x = y (m/pixel) that provides the DTM scale and temporal discretisation according to a time step t that ensures the stability of the computational process, verifying the Courant-Friedrich-Lewy (CFL) condition [131] ||u||∆t/∆x << 1 [132]. The CFL condition, necessary but not sufficient for the stability of the numerical scheme, states that the distance travelled by any mass at velocity u during the length of the time step ∆t within the mesh ∆x must be smaller than the distance between the mesh elements.

CFL Adjustment Condition for Digital Terrain Model
In the simulation process, the resolution and content of the topographic information incorporated in the calculation through the DTM play key roles in the evolution of the plan-view distribution of the lava site [33,133,134]. Generally, low-resolution DTMs provide poor results that do not match the actual perimeters of the lava flow site [47], especially in areas where the natural topography has been highly modified by humans [135]. Additionally, terrain roughness affects mobility as terrain slope variations intervene in the above SWE system. So, to assess the possibility of constructions or other roughness in models on the lava flow simulation [4,136] in real time, in this study, we used used two types of DTM, already discussed in the Materials section: (i) DEM containing the terrain elevation data, and (ii) DSM containing relevant buildings without vegetation, taking into account that all buildings equally behaveas terrain protrusions without considering their constructive characteristics.
Therefore, before carrying out the calculation in this study, the DTM resolution was selected at ∆x = ∆y = 6 m. From this value, the time step size was set so that: (i) it was small enough to verify the CFL stability conditions, and (ii) it was high enough for the calculation to be performed as quickly as possible on a computer. With an appropriate value, we ensured that the mass of a cell or grid element of the spatial discretisation, set by the size of the DTM, propagates at most at velocity u only to the immediately neighbouring cells.

Calibration Criteria
The purpose of the calibration of the planned scenarios to simulate ( Table 2) was to estimate the set of rheological parameters that make it possible to obtain results similar to that observed on a DTM and for a specific period of time. The diversity of scenarios examined in this study (SS1 to SS8, Table 2) aimed to capture the possible factors, other than the physics of lava movement, that can affect the geometry of the lava flow emplacement. One of such factors is the morphology of the terrain and its resolution. Therefore, we simulated different scenarios including and excluding buildings and other obstacles (trees and bridges). The DSM scenarios (SS4, SS6, and SS8) included obstacles, whereas the DEM scenarios (SS1, SS2, SS5, and SS7) used bare terrain. Two other factors are the location of the volcano vent and the amount of lava sequentially emitted by each vent. So, in case there are several sequential vents, the simulation was performed in several phases, as many as there were vents, with each one simulating the outflow from each vent over the previously emitted lava.
The similarity between a simulation and reality is usually qualitatively and quantitatively evaluated on the geometry of the lava flow and, when data are available, also on the thickness and velocity of the lava flow. The qualitative option is only suitable for monitoring large-scale behaviour, comparing the order of magnitude of reaches or the extent occupied. For quantitative evaluation, a posteriori metrics are used, such as those proposed by [67,133,137], among others, which can be calculated at the precise time steps of the simulation matching the time for which ground truth geometric information is available.
Parameters can be semiautomatically estimated by numerically solving an optimisation problem that poses the minimisation of an objective function based on the distance metrics between the real and simulated lava flow used for quantitative comparison [138]. The choice of an appropriate metric function is essential to evaluating the goodness of fit of a simulation. Due to the number of parameters to be calibrated and the degree of nonlinearity of the problem close to the minimum, the computational cost of finding a best fit can be high. The option of manually determining an approximation of the minimum was considered valid and fast enough for this study.
To quantitatively verify the correspondence between observed and simulated lava flow lengths, some researchers [78,[138][139][140] have used the length rate S L as a metric, which is defined as: where L A and L B are the real and simulated distances, respectively, taken from the vent point to the front of the lava flow in the direction of the main lava stream (Figure 3). S L quantitatively evaluates the performance of the simulations by comparing the simulated with the observed run-out using the length ratio. Likewise, to express the degree of similarity between the real and simulated lava extents, we used the Jaccard coefficient [141][142][143] when given two nonempty finite sets A and B: which heuristically measuresthe probability that an element of at least one of the two sets is an element of both. Thus, S O is a reasonable measure of the similarity or "overlap" between the two (Figure 3). The numerator in (14) is the area occupied by both observed and simulated flows and is the true positive. In the denominator of (14), the area filled only by the simulated flows is the false positive or overpredicted area. The area covered only by the observed flows is the false negative or underpredicted area. Therefore, a value of S O = 1 is a perfect match between actual events and simulations. Although some versions of this coefficient exist [138], the Boolean similarity measure S O (16) has been frequently used in other studies related to the comparison [67] calibration of lava flow emplacement models [70,134,[144][145][146], including to compare the effects of DEM resolution [133] or the differences of simulated results on a DEM or DSM [135].
In this study, we set A and B in (13) and (14) referring to the surface occupied in reality and the simulation resulting from the scenario model, respectively, both taken on the same date ( Figure 3). Therefore, the values of the rheological parameters were manually modified, controlling the behaviour of S O throughout the days of the simulation (from 19 to 25 September), trying to achieve the highest value of S O at the end of the simulation. This iterative procedure was applied for both the simulated scenarios ( Table 2) with one and two sequential emission vents.

Simulations Set-Up
To run each simulation of the lava flow in the period of the La Palma eruption considered in this study (Table 2), we incorporated the estimates of the rheological parameters (ρ, η, ϕ int , ϕ bed , R 0 , ξ) obtained from the analysis of the lava flows mapped in the field by COPERNICUS (Table 1) and from the bibliographic sources of the VolcFlow software authors (mainly R 0 , ξ). It is assumed that, for the short period of time analysed, the lava emission rate Q was constant and was obtained by evaluating the volume of the lava flow over time since the eruption started. The topography of the pre-eruption terrain was incorporated in the calculation by means of a raster DEM elevation or DSM surface file with a sufficiently high resolution and on which the effusive sources were marked.
The  (Figure 2). As the sources changed during the study period, the emission flux was distributed between each of them proportional to the time they were active. A precise map of lava vents is not required provided lava flow concentrates into a single flow tongue very close to the vents. Hence, providing the overall emission rate suffices to appropriately consider the flow. In addition, the DTM on which the simulation was performed was modified according to the height distribution resulting from the location of the first source. Therefore, the simulations including the second sources were carried out on modified terrain according to the emplacement of the first lavas. This second source was placed following the same assumption as the first one, not requiring a precise mapping of its location rather than appropriately capturing its emission rate toward generating the flow.
Along the complete processes, between one source and another, the geometry of the lava flow was evaluated against the mapped geometry following the calibration criteria to readjust the rheological parameters of the lava. With this first readjustment, the calculation was relaunched from the beginning to the end of the simulation time. The result was re-evaluated with the information on day 25, when the thickness distribution was also available. Numerically, to maintain the necessary condition of CFL stability during the entire simulation, and as the DTMs used here (DEM and DSM) were taken with ∆x = ∆y = 6 m, an integration time step of ∆t = 0.5 s was adopted, according to the maximum velocity observed throughout the evolution, a time step approximately two times lower than needed by the CFL condition. Below this resolution in the DTM (although the original data allowed better than 1 m resolution), computing time grew dramatically because, in addition to increasing the number of points over which the numerical approximation of (7) to (9) was solved, to maintain the CFL condition, it was necessary to take a much smaller ∆t. Above this resolution, even though the computation time decreased, the representation of the relief loses obstacles and roughness that were averaged or smoothed, being critical in the displacement of the lava flow, were not properly represented.

Flow Rate and Rheology Estimation from Morphometry
Information about lava thicknesses evaluated on 25 September, as a result of the optical UAV flight carried out by the Cabildo Insular de La Palma (Table 1), allowed us to evaluate some descriptive thickness statistics, such as its sample distribution or cumulative frequency function (CDF, black line in Figure 4), its median h 50 = 13.65 m, its average <h> = 12.04 m, and standard deviation 9.5 m, with values up to 30 m at some accumulation points. h 50 > h is indicative of the asymmetry and positive skewness of the distribution function, meaning that higher thickness values are less probable than lower values. This skewness in the distribution requires that estimates of the lava volume and lava viscosity, using Equations (1)-(3), and emission rates Q, using Equations (4)-(6), consider the possible range of variation due to skewness. Therefore, when calculations were performed with these equations, the value of h was replaced with the median h 50 or the mean <h>.
The morphometrically estimated averaged and instantaneous velocities showed an overall fairly constant behaviour, with a slight decreasing trend, over the first 5 days since the onset of the eruption (dashed curves in Figure 5a). The average values obtained (filled symbols in Figure 5a   Graphically (Figure 5b), the result of applying the different approaches introduced in the Method sections to the calculation of the volume of lava being emplaced in the study area showed that, in all cases, the volume showed an increasing monotonic behaviour. At the beginning, the amount of emplaced volume more quickly increased than at the end of the time period analysed, as shown in simplified form (Figure 5b) in the dashed curve. The volumes estimated using the [WL] i data (Table 1) from the COPERNICUS (void symbols) and the Cabildo de La Palma (full symbols) are very similar, with a value close to 3 × 10 7 m 3 after the seventh day. The differences in the results, whether the mean or median heights, that were used in the calculation were not really significant. However, if the frequency distribution of heights was used, although the behaviour over time was very similar, the estimated values were lower than the former, reaching 2 × 10 7 m 3 , i.e., estimating one million cubic metres less.
The flow rates Q obtained from the previous results of the volumes emitted and the descriptive statistics of the lava heights, by applying Equations (4)-(6), showed a behaviour similar to the eruptive process evolution (Figure 5c), with an increase in the first two days, followed by a decrease and then a tendency to be constant over time. The maximum rates estimated by expression (7) using h 50 was of the order of 140 m 3 /s (1.42 × 10 2 to 1.36 × 10 2 m 3 /s) that, from 23 September, began to decline toward 55 m 3 /s, using expression (5), to 43 m 3 /s, obtained with the average <h>, and expression (4), with small variations in time within this interval around 50 m 3 /s on average, which were considerations apart from the source of origin of the data used, which fell within this oscillation. However, the Q values obtained for the last days considered in this study, using the volume estimated with the histogram of the height frequencies ( Figure 4 as the cumulative probability distribution function) and expressions (4) and (6), were about 30 m 3 /s, lower than the first ones estimated with <h> and h 50 .
Finally, the results obtained on the rheology of lava, through the viscosity η (Figure 5d), showed some disparity with each other depending on whether Equation (1), (2), or (3) was used. Using Equation (1), we obtained (points in Figure 5d) an increasing trend behaviour for the viscosity, even when using the mean or median of the lava heights (as shown by its dashed curve), whose values ranged from nearly 7 × 10 6 Pa s up to 3 × 10 7 Pa s at the end of the period observed in this work. Despite the different behaviour of Equations (1) and (2), both seem to converge over time toward a similar viscosity value in the range of 3-4 × 10 7 Pa s, although it would require more observation time to determine whether they remain the same or diverge from this point onward. Higher values than those plotted in Figure 5d would be obtained if the emission rate obtained from the frequency histogram of the heights was considered, because Q appears in the denominator of (2) and (3). However, the viscosity response obtained with Equation (1) increases as function of time and independent because it does not play a role in it.

Numerical Simulation and Calibration Results
The values of the emission rates and viscosities that we obtained morphometrically were used as initial parameters in the numerical simulations (Table 3) for each scenario. They were also used as initial values in the iterative optimisation process, minimizing the Jaccard coefficient through manual adjustment, and calibrated for a more precise estimation. According to the proposed simulation plan (Table 2), we considered the case that the SS* scenario was performed in a single phase, i.e., only one vent (just V1 in SS1, SS3, SS4, and SS8) throughout the simulated time: the asymptotic average flow rate Q value (50 m 3 /s, Figure 5c) was assigned as the initial Q1 value for this point for the optimization process. If two or three simulation phases were considered in the SS* scenario, that is, considering two vents (V1 and V2 in SS2, SS5, and SS6) or three vents (V1, V2, and V3, in SS7) vents, depending on their reported onset times (Table 2), the average Q for the first and second days (131 m 3 /s, Figure 5c) was taken through the first vent, and the asymptotic value (50 m 3 /s) for the rest of the simulated days, depending on the time scope of the scenario to be simulated, was taken through the second or third vent. In scenarios with two or more phases, switching from one phase to another, the lava source point vent changed and updated and modified the DTM with the lava emplaced up to the time of the phase change.
Concerning the viscosities to be initially incorporated in the optimization, the same strategy as above was followed. Assuming that within the simulated time span, the viscosity increased due to heat loss, if the scenario to be simulated involved a single vent, the initial viscosity value was 3 × 10 6 Pa s. For scenarios involving two or three phases, the viscosities initially used assumed the estimated values for each phase: 7.6 × 10 6 Pa s for the first phase, and the estimated asymptotic value of 3 × 10 7 Pa s for the second or third one. The rest of the rheological values (expression (14)) used to carry out the simulation of the numerical model (ϕ int , ϕ bed , R 0 , ξ) were obtained from the values recommended in the literature [33,126,128] that made use of the VolcFlow program for the simulation of lava flows (Table 3). Then, an initial internal friction angle (ϕ int ) limited by the average slope of the terrain between the point of emission and that reached on 27 September (10 • ), a density of 2100 kg/m 3 , and a constant Voellmy coefficient ξ = 0.01 were taken for all scenarios. Consequently, in the results obtained from the calibration exercise for each scenario, we present the best-fit approximations for Q and η: the values of the internal and basal friction angles and yield strength R 0 that provide the best values of the Jaccard coefficient. Table 3. Arbitrary initial lava emission flow rate Q (m 3 /s) from vents, and apparent viscosity η (Pa s) from morphometric analyses used for calibration through simulation. Simulation flow rate Q and rheological characteristics (viscosity, friction angles ( • ), and yield strength (Pa)) estimated from calibrations in each scenario (SS* in Table 2). The simulations of the first five days performed in each scenario (SS1 to SS6, Table 3) were calibrated against a Jaccard coefficient that, in the best case, was higher than 70% in SS3 and SS4 ( Figure 6). The rest of the scenarios were calibrated with Jaccard coefficient values greater than 60%. In all these cases, and regardless of the digital model on which the simulation was performed (DEM or DSM) or the Jaccard coefficient value achieved, the lava height distributions obtained by numerical simulation on 25 September were very similar to each other and to the real one on the same day for values higher than 15 m (Figure 4). However, the values of lava thicknesses below 15 m did not sufficiently match the real distribution, with the difference being greater as the thickness was smaller, as shown in the cumulative probability function (Figure 4), to be comparable with the available field data as originally distributed by the Cabildo de La Palma.

Simulated
The resulting fit between the simulated flow fields with VolcFlow in each scenario and the actual observed flow field was very close, but not exact. It appears that the flood extent was well-bounded laterally, but underestimated at the advancing front, regardless of whether one or two phases was used in the simulation. These differences in space and time were probably caused by the relatively high uncertainty in the initial morphometric estimates of the emitted flow rate Q, the viscosity, and the nonuniform distribution of the flow along the eruptive fissures. Also critical for these simulations was the unavailability of DEMs that are frequently updated to account for topographic changes caused by ongoing eruptive activity (both lava flows and cinder cones) and of a resolution high enough to simulate the expected activity. As expected, the scenarios calibrated with a low-resolution DEM (SS1 and SS2) were the ones that produced the lowest Jaccard coefficient values, about 63% to 65%, being somewhat better (at the beginning) if two effusive vents were used than if one was used. However, later, after 22 September (Figure 6), they were very similar and did not have the expected length (Figure 7a,b). Specifically, despite the differences in the maximum lava flow distance, the effect of using a 6 m/pixel DSM (SS4, Figure 7d) significantly improved the result, obtaining a Jaccard coefficient of 72.5%, compared with using a DEM with the same resolution (SS3, Figure 7c) that maintained a difference of 4% over the course of the simulation period. The effect of introducing two vents, trying to reproduce the sequence of the observed process, did not result in an improvement in the fits to the real lava extent (pink areas in Figure 7). In these cases (Figure 7e,f), the uncertainty of not knowing the part of the flow that flowed from the second one, to be distributed between both, although resulting in a relatively high Jaccard coefficient (68% to 69%) for either DEM or DSM, the range was not sufficient, which is key for the simulation to provide useful information for early warning.
The best fit provided (last three columns in Table 3) by the numerical simulations preformed for each scenario up to 25 September allowed us to estimate the values of the flow rate and the rheological parameters of the lava according to the assumption of one or two simulation phases. For a single phase (for 20-25 September; scenarios SS1, SS3, and SS4), the best-fit flow rate was 65 m 3 /s. However, if two phases were considered (scenarios SS2, SS5, and SS6), the rates considered the initial intense transient eruptive process and the subsequent descent into a sustained regime, so that, for the first 24 h of emission, 140 m 3 /s was the best estimated fit. For the remaining 4 days, 58 m 3 /s was estimated. In the numerically calibrated viscosities, a different behaviour was also obtained, 2.9 × 10 5 Pa s if a five-day phase was considered, or 2.4 × 10 3 Pa s for the first 24 h and 2.9 × 10 5 Pa s if the calibrations were made with simulations in two phases. The rest of the rheological parameters obtained: yield strength and internal and basal friction angles, showed no differences in the adjusted values (Table 3) regardless of whether they were simulated in one, two, or three phases.
The calibrated effusion rates and rheological parameters were used to generate a longer-term evolution forecast until 27 September ( Figure 8). As expected, despite the results previously obtained on the DEM and incorporating all the sources that arose in the field during the simulation time (scenario SS7), the forecast result obtained for 27 September (Figure 8a) did not achieve the observed range. Although it laterally exceeded the limits of the perimeter of the real lava flow (area in green), its Jaccard coefficient was far (57%) from an acceptable target. Finally, the SS8 scenario presented a good fit between the model output and field data (Figure 8b). The numerical simulation also included the contribution of one of the ephemeral vents: the first detected and located (Table 2), with satisfactory agreement between the model output and field data ( Figure 6), which reached a Jaccard coefficient of more than 85% over a DSM of 6 m/pixel and a distribution of lava flow heights similar to the real ones recorded in the field for the highest values.   The time required to undertake the simulations includes data gathering, data processing, and CPU time to reach a numerical solution to the model through the entire calibration process; hence, it is the actual time spent working on the model. The simulation time interval is the range between the date of the year when the process to be reproduced by the model starts and ends. Thus, the CPU time in Figure 9 represents the amount of time the computer required to perform a calibrated simulation in the simulation time interval. Therefore, anything below the 1:1 line means faster than real-time simulation was achieved. The times in Figure 9 were obtained using VolcFlow running on MATLAB ® R2021a (Math-Works, Inc., Natick, MA, USA) using a computer with an Intel ® CoreTM i9-7900X CPU @ 3.30 GHz with 128 Gb of RAM (Intel Corp, Santa Clara, CA, USA). Computation time increased with the resolution of the digital terrain model: we obtained 0.28 days with a 12 m pixel DEM up to 2 days on a 6 m DEM. The complexity of the digital terrain model also influenced the calculation time, reaching up to four days using a 6 m DSM (terrain with obstacles such as buildings). The ratio of CPU time and the time span simulated improved as the simulation span increased. A 24 h forecast was achieved for a simulation span of five days with a DSM with 6 m resolution, and a 48 h forecast was achieved with a simulation span of seven days. These results are a promising step toward real-time forecasting as it can already be achieved using the calibration in this study for a similar scenario. It also means that manual or automatic calibration would require much faster computation depending on the number of runs required to undertake such calibration, with a rough estimate of at least 10×. This study did not include an analysis of the converging results of the successive iterations assessed manually here.

Discussion
Emergency control systems face numerous challenges during volcanic eruptions. Preparing hazard maps [49,66,79,144], monitoring the progress of the lava flows [33,40,70] and forecasting lava displacements are amongst such challenges [47,137]. For this purpose, the use of numerical models is a great development, particularly because they can be run increasingly faster using more efficient algorithms and more powerful computers [86,148]. When using physically based numerical models, it is necessary to take into account the set of equations and simplifications they solve, which are derived from the equations of conservation of mass, momentum, and energy [67,68,128,129,140]. Regardless if the strategy follows cellular automata [70,72,138], finite differences, elements, or volumes, [77,149,150], they all need input data that are usually not easy to obtain [61,103,115]. This difficulty may be prevent the provision of forecasts based on models. Our study proves we are on the verge of a change in paradigm, moving from expert assessment and using previous maps towards providing forecasts. This paradigm swift requires specific measurements and observations during events and their immediate public release.
Our study is not exempt of the uncertainties derived both on the simplification of the conceptual model of lava flow and the simplifications assumed. The optimization criteria (the Jaccard coefficient) may also have contributed to the uncertainties, as other optimization criteria could have been used. Morphometric measurements of the lava flow were also subject to uncertainties due to spatial variability and measurement errors owing to the resolution of the image used, the calibration of the measuring equipment, and possible accidental errors that may have randomly predisturbed the data. The information available in open access and that has been used in this study to perform the morphometric analysis did not allow us to evaluate the uncertainties in the width, length, height, and slope of the lava flow and its propagation in the estimation of volumes, emission rates, and rheological parameters, as in [121] for small portions of the flow. The variance in the thicknesses, in addition to being a statistic referring to the whole lava field, is a value that refers to the variability in all the thickness values and does not estimate the uncertainty of each one. However, the values of volumes, emission rates, velocities, and thickness of lava were provided with sufficient accuracy; and the time to provide forecast is already within reach of actual technology. These results are not only numerical, but also morphometrical. Regarding the morphology of lava, sometimes our approximation overestimated the observations, whilst on other occasions, our results were underestimated. Overestimation occurred when approximations followed instantaneous velocities and average thickness. Underestimations may have responded to the conceptual model, where very fluid lava flows resulted in a set of anastomosing lava tongues and channels, not to mention lava tubes or younger flows overlapping older, yet still active, flows. Observations accounted for the maximum accumulation thickness of about 70 m and average values of 12 m throughout the eruption [51]. Our model provided very similar average values, including that of individual lava tongue thickness, due to the constraints imposed on the model and its calibration.
In the initial stage of the eruption, our results showed that emission rates, if obtained based on averages (Equation (4)), were around one order of magnitude larger than those obtained using average velocities of the front and weighting according to the histogram frequencies of thickness (Equations (5) and (6)). Although these initial differences decreased as the eruptive process evolved, on the sixth day, for the Tajogaite event, the estimates were very similar. These initial differences may have been due to the fact that the velocity in the lava flow was in a very disparate transient state, and the velocity field measurements were not able to capture the average behaviour, producing underestimated values.
Recently, [51] estimated a 27 m 3 /s the emission rate and 50 × 10 6 m 3 [11] of erupted material within the first eight days of the eruption (21-27 September), while we estimated an average asymptotic emission rate of about 57 to 60 m 3 /s and, 30-35 × 10 6 m 3 , respectively. Both values are well aligned considering the lower reflects the mean effusion rate of the eruption [88,115], while the larger values represent the average during the first five days of eruption. During the first days, the event was still in the transience of the beginning of the eruption, without having reached a phase with a stationary flow rate. Despite possible improvements in the estimation values, they should be checked with other estimation tools independent of morphometry [115]. An interesting one is based on the maximum estimate of the time-averaged discharge rate (TADR), calculated [151,152] from satellite data in the infrared band (MODIS, AVHRR, SEVIRI, etc.). However, converting from satellite thermal data to TADR requires the use a series of physical parameters of the lava (density, heat capacity, vesicularity, latent heat, etc.), which are not easy to obtain in situ during the eruption and must be taken from literature sources [115]. Therefore, the validity of this remote estimates is also uncertain. Comparison of satellite-derived discharge rate estimates with those derived morphometrically will allow the derivation of appropriate values and their variation over time as an input to the numerical simulation models [40,90].
The lava viscosity obtained in this study using Equations (1) [110] and (2) [111], assuming the Newtonian behaviour of the lava, incorporates neither information on the composition nor on the crystal content of the lava nor its temperature. Similar to what occurred in [121], the results of these two morphometric approximations are admissible proxies in a stationary regime. They approximate viscosity at a particular time accord-ing to the observed geometry. The evolution of viscosity behaved in a decreasing way when Equations (2) and (3) were used, being slightly faster with Equation (3) than with Equation (2) (see dashed curves plotted for each result in Figure 5d). The results were acceptable when the lava flow did not develop lateral levees and its extension was considerably greater than its thickness [117]. This was the case of the first days we studied. However, the morphometric information used here was taken almost in real time during the eruption, with a one-day lag (the time COPERNICUS needed to provide their maps), and not once the eruption had finished, as used by [121]. On the other hand, the values given by expression (3) significantly differ from the two previous ones, as they did not incorporate the effect of the slope or the lateral extent of the flow, which were important in this case. However, the results of the viscosities obtained with formulae (1) and (2) differ by one to two orders of magnitude from those obtained by numerical calibration for the asymptotic rheological behaviour of lava. The Jeffreys equation assumes a Newtonian behaviour, but if significant cooling has occurred, the viscosity also increases. Recent studies in the area to analyse the viscosity of lavas from the Cumbre Vieja eruption have shown them to be an extremely fluid and fast nepheline normative basanite [51,65], alkaline with low silica (<45% SiO 2 ), which produced cascading lava flows and a remarkable development of an expansive lava flow field that took only a few weeks to form. These phenomena are distinctive of a rheology with viscosities exhibiting values from exceptionally less than 1 × 10 0.83 to approximately 1 × 10 5 Pa s (rheometer-based measurements) within eruption temperatures of about 1300 to 1100 • C [65], respectively, with a step change to higher viscosities at T < 1125 • C, highlighting the influence of high crystal contents. Compared with those obtained by our numerical calibration (Table 3), they are an order of magnitude lower Although [65] these are exceptionally low viscosity values for high temperatures, they are not far from the values due to low silica content [103]. This would be particularly true if we take into account the important uncertainties associated with morphometric estimations and manual calibrations. Furthermore, if the viscosity rheometer measurements refer to point data in lava rivers and not the average behaviour of the lava as a whole, a continuous media would be numerically modelled through Equations (7) to (9), whose velocities are much lower (10 −2 m/s) than those measured in lava rivers (10 m/s according to [65]), which presented cascades. We would observe abrupt transitions of a flow from supercritical-to subcritical-flow (hydraulic jumps) states, or high-speed channels with Reynolds numbers close to turbulent regimes, not provided with the system of partial differential Equations (7)-(9) used in this study.
The viscosity values estimated here from the morphological dimensions of the flow are quite low compared with those obtained in andesitic lava flows (10 11 to 10 13 Pa s), and do not significantly differ from the values estimated from field data. The value obtained is slightly higher, as it is a value averaged over the complete extension of the lava, which is stopped by the important surface cooling with an important crust. However, in its interior, there are still important molten material flows as a result of its low viscosity. The content of crystals in the whole of the lava corresponds to the rheology of the flow, in which it becomes governed by the solid particles and marked by a slight change in slope, as shown in Formulas (1) and (2). In its displacement, the lava reaches the rheological limit, where it is too viscous to continue flowing (for this given strain rate). This suggests that the flow, as a whole, has registered the highest viscosity when it is apparently stagnating in its morphology. Therefore, as [121] suggested, the rheological values estimated from the morphology are very sensitive to the crystalline fraction and, above all, to the crystal shape (and, hence, the crystalline phases). In the context of the Cumbre Vieja basanitic flow, as it cools and crystallises, the viscosity increases until it reaches the value corresponding to the time scale for which the flow appears to be stagnant. In other words, the morphometrically derived viscosity of a lava flow is a promoted and instantaneous apparent value that is defined by the morphological dimensions of the lava flow, resulting from the irregular displacements according to the nonuniform coupled conditions of temperature, crystallisation, viscosity, and yield strength that occur at the different advancing fronts.
We validated the results obtained with VolcFlow code by comparing some morphological characteristics of real-world flows (e.g., length, area, and thickness) against those generated by the simulations [76,79,97,128]. These comparisons helped to improve the results of numerical calibration. They also aided in the understanding of the complex rheological, thermal, and dynamic processes involved in lava flow emplacement and, therefore, in more accuratelysimulating lava flow emplacement events [129]. They also shed light on whether the given assumptions still allow the reproductino of the actual emplacement of the flow. In particular, the important differences that arise when performing simulations on a DEM and on a DSM helped to identify the effects of model resolution and artificial structures on the evolution of the lava tongue [4,133,135]. Very locally, on small tongues, it may be necessary to detail the properties of the thermal behaviour of these structures against the advancing of the lava. On other cases, infrastructure, such as greenhouses incorporated into the DSM, acts as a true barrier even though their resistance to lava is minimal.
Because ϕ bed is affected by several parameters such as the basal surface roughness, flow thickness, grain-size distribution, topography, channel irregularity, presence of obstacles or hydraulic jumps, among others, and additionally, because most of these parameters strongly vary along the flow path, it is difficult to find a domain-wide ϕ bed capable of capturing the complex interaction of these parameters and their local variations. Its assessment becomes even more complicated when preliminary information on the past deposit is not available. As Equation (10) was derived assuming simultaneous Mohr-Coulomb behaviour at the contact of the lava flow with the ground (ϕ = ϕ bed ) and within the lava (ϕ = ϕ int ), when ϕ bed = 0, Equation (10) reduces to the classical Rankine definition. The results fitted with the model suggested that the lava mass is mobilized with a rheology that combines the Mohr-Coulomb lateral stress ratio with the Bingham behaviour. The zero basal friction angle implies that the basal Coulomb shear stresses are negligible (R coulomb = 0) with respect to the other shear stresses involved in R (Equations (11) and (12)). From the Mohr-Coulomb perspective, once the lava has been emitted, it moves along the slope when the slope is greater than the basal angle, that is, from the first moment and by the effect of gravity (first term of the second member in Equations (8) and (9)), without the retarding effect of the terrain roughness. As the value of the internal friction angle is also very low, close to the basal value, it indicates that the movement of the mass is almost in the form of a slab of variable thickness. Whether the basal friction angle has a negligible effect on the motion, it may indicate that the lava flow behaves, in addition to the Rankine behaviour due to the nonzero internal friction angle, as a Bingham fluid, i.e., a combination of the rheological models of a viscous solid and a plastic fluid, with added behaviour representing the effect of turbulence and/or collisions as a function of the lava flow velocity [126]. The plastic rheology of lava was characterised with a yield stress limit R 0 of between 35 × 10 3 and 42 × 10 3 Pa, a value independent of velocity and thickness that limits mobility, because when the drag stress drops below R 0 , the material decelerates and stops.
Despite these drawbacks, the simulations carried out in this study, for the seven days of advance over the DSM, very accurately reproduced the geometry of the lava flow (85% Jaccard coefficient). The best-fit morphology was that of greater thicknesses, which, in the end, can cause more damage to infrastructure. This also allowed us to calibrate the average values in the rheology of the fluid, which were very similar to those estimated by other means and just recently published. The calculation timing, despite the complexity of the problem, and without being computationally improved, is adequate to make a preliminary estimate based on morphometry or other experimental rheometric in situ techniques [65,103], and to improve this first estimate by semiautomatic numerical calibration. This calibration may provide a better and complete rheological characterisation and a short-term forecast of the evolution of the front.
Our study, comparing the performance of the model using bare terrain (DEM) versus terrain with obstacles (DSM), provided results similar to those of other authors. We also achieved better Jaccard coefficients for DSMs, but the difference we found was not as notably high as that of [135]. Some obstacles behave like dams (road talus), others like channels (road clearances), and others like small pits (water rafts). However, not all infrastructures should be integrated into the DSM for lava flow modelling purposes. Greenhouses (very abundant in our test site) and light walls show very low resistance to lava fronts, but the DSM accounts for them as true barriers that locally disturb the flow. Our study involved lavas flowing over towns; yet, flow obstacles were not very close to each other, as opposed to of [135], which, in our model, led to softening the effects of the obstacles. At the same time, we also found that the higher the resolution of the DTM, the better Jaccard coefficient. Nonetheless, using very-high-resolution DTM had a deep impact on model performance. The CFL condition would not be met even with very low ∆t, and the time required to increase model instability was already too long to be considered for forecasting purposes. A ∆(x, y) of 6 m in the DSM was capable of keeping most of the relevant infrastructure, as well as erasing most of the smaller items, such as walls or sheds. With lower-resolution models, the impacts some items decreased, such as the many roads traversing the area. Still, this effect did not considerably impact the goodness of fit.

Conclusions
Our study should inspire the community to research lava flow forecasting, as we proved it may already be possible with our current technological advances and data gathering capabilities. Moreover, our study contributes to further understanding lava flows and may encourage an increase in model complexity hereafter both for forecast purposes as well as for planning and mapping lava hazards.
We calibrated and validated our model with limited real-time data. Our average lava flow rheological fits, despite the simplifications of the numerical model used in the calibration, within the range of the experimental evidence for the measured viscosity of the alkaline lava based on the Tajogaite eruption. The observations of the morphometry distributed in time is almost as simple as possible for validation and calibration. Including lava thickness was a beneficial addition, and we call on the monitoring community to provide open-access data that also account for lava thickness. The differences imply that the geometric approach of viscosity cannot accurately constrain the details of the rheological evolution of a flow, although it may be a valid estimate.
The results obtained in this study can certainly be improved, for example, by including less simplifications or using more complex approaches, such as those presented in the Introduction. Other codes should also be included in similar experiments to further test the robustness of the results. With or without improvements, during a volcanic crisis, a balance must be found between precision and forecast time. In Spain, the civil protection authorities have repeatedly suggested that early warnings may be effective if dispatched at least one hour in advance of the hazard. Due to the nature of lava flows, our findings suggest that lava forecast will be feasible very soon. Therefore, in case of an emergency, a quick response time is more essential than the delay that can be caused by working with complex models that demand a large volume of data for their calculation. With the minimum amount of data available for an isothermal model such as the one used in this study, it was possible to sequentially incorporate the changes in the volume emitted and in the rheology that would be produced by the cooling of the lava. The use of a model with cooling [128,153] implies a higher computational cost, having to solve the energy balance equation at each time point, and a greater number of input data need to be adjusted by optimisation, as they include the variations in viscosity and yield strength with temperature. In this case, in order to reduce uncertainties, the optimisation problem may be formulated with an objective function (or fitness function) that not only minimises the geometric Jaccard distance, but also takes into account other metrics between the real reference casting (described as a polygon) and the simulated one, such as those used in [137], and between thickness distributions (lava flow height variation within the reference polygon) and flow velocities. This information can be collected in the field, by satellite, or by UAVs.
As the mathematical formulation of this problem with a multiobjective function (geometry, thicknesses, and velocities) is hard to numerically solve, a semiautomatic calibration through hundreds of simulations is required to locate the minimum of the problem. Therefore, to answer the question posed, we cannot claim that it is actually possible to achieve real-time lava flow estimation of sufficient quality to inform a volcanic emergency response, although we are on the way to further improvements. By real-time estimation, we mean the workflow that obtains the lava inundation area at a future time. This workflow, when the behaviour of the lava is not known because it is a new volcanic cone (as we presumed in this study), consists of the following: (1) the initial estimation of a rheology from the first instants of lava emplacement, which can be achieved using the morphometric methods applied here or those suggested in [121]; (2) the fit of the rheology by optimisation using geometrical (perimeter and thickness) and kinematic (flow velocities) constraints with an appropriate numerical model to capture the thermo-hydraulic behaviour of viscous-plastic lava in its equations; and (3) the simulation, using the properties obtained in the previous step (2) of the forecast of the lava extent at a future time. Although recent studies evidently made efforts to achieve real-time simulations, these only referred to the calculation speed of phase (3) being higher than that of the lava flow advance in order to achieve on-time future predictions. All of them, being performed on well-known cases, assumed that the lava model input properties are previously well-known and well-determined in phase (3). In particular, in [66,[153][154][155], the emission rate from the volcanic cone that feeds the lava flow was estimated, a time function that constituted the source term of the equations. This function was considered to be the only particular feature of the ongoing eruption that was estimated from remote satellite data and was used as a real-time input to the numerical lava flow model. On the other hand, [153] introduced a depth-averaged finite-volume model, considering the rheological volume changes due to thermo-rheological stratification, allowed the estimation of more realistic vertical lava profile rehology. From the perspective of our study, this recent advanced model may have interesting applications in analysing and reproducing the complexity of the lava flow in a posteriori scenarios. Particularly when the large set of physical parameters required for its use is known. However, its usefulness for forecast operation in active volcanic regions might be limited due to its computational efficiency and the intensive demand of data.
Reducing calculation times may be possible by improving the computational performance of the software, working in a parallel computer environment, and vectorising their processing so that they are carried out mostly on GPUs. On the other hand, the institutions in charge of tracking and monitoring need to know the information that the models require, so they may be able to sample and collect these data. Providing these data in open access will greatly contribute to the scientific community and decision making.
We encourage the community to use full digital terrain models (DTMs) that include buildings and infrastructure (DSM). They produced better overall performance than bare terrain (DEM), even if low-resisting infrastructures were present in the DSMs. We also suggest lightly preprocessing the DSMs when performing lava modeling to erase as many of those softer obstacles to lava flows as possible, such as greenhouses or even light industrial units. If modelling time is not an issue, it seems that the higher the resolution, the better results. There is a limit to this assumption, as the more detailed the boundary conditions, the more complex the models. Still, a balance can be easily found between computing performance and goodness of fit. A DSM of about 6 m resolution produces good balance; however, for automatic calibration, computing intensity may be too high without intensive parallelization of the code or better computing resources than those used in our study. Running the code in GPUs, as it has been performed for other types of geophysical fluids, should also be explored. Funding: This work has been partially founded by Project COGNOSCI-ON: Generation and assimilation of knowledge (COGNOS) through scientific activation (SCI-ON), in subjects related to engineering in Earth Sciences, of the Competitive Program of Educational Innovation at the Universidad Politécnica de Madrid PIE-UPM (IE1920.0603) through 19-20 and 20-21 grants. This work has also been supported by the "Severo Ochoa Extraordinary Mention Funds" from IGME-CSIC.

Acknowledgments:
We are grateful to Karim Kelfoun for his kind help in providing support with VolcFlow software, and to all the reviewers whose comments have greatly improved the manuscript, guiding its development to its final version. We also want to thank all public authorities that provided data in an open manner and, above all, we would like to thank the emergency responders and collaborators for bringing our families and friends in La Palma to safety. We are deeply moved by the resilience of those who are still awaiting solutions.

Conflicts of Interest:
The authors declare that they have no known competing financial interest or personal relationships that could have appeared to influence the work reported in this paper.

Abbreviations
The following abbreviations are used in this manuscript: