A lava flow simulation experience oriented to disaster risk reduction, early warning systems and response during the 2021 volcanic eruption in Cumbre Vieja, La Palma

Lava flows are one of the hazards involved in a volcanic eruption, and although they rarely cause the loss of human life, they are highly destructive in terms of damage to property and economic activity. Therefore, the management of volcanic disasters requires fast and accurate information on the behaviour and evolution of the flows, mainly related to their extension, displacement, and trajectory. This was the case during the disaster linked to the volcanic event that occurred on the island of La Palma in the Cumbre Vieja area at the end of 2021, which lasted eighty-five days. This paper describes part of the work performed by many different groups to provide predictive information aimed at feeding the early warning system set up during the disaster. This case shows the experience in the use of a probabilistic simulation algorithm implemented in the Q-LavHA plugin for the QGIS software, which is both easily accessible and applicable, to analyze its features in detail, as well as its predictive capacity. The results show that the model can efficiently and quickly satisfy the demand for this type of information, and its high similarity value is also validated by the Kappa index.


Introduction
Lava flows are one of the main hazards during the course of effusive volcanic eruptions posing a risk to the population, and therefore require constant monitoring. In general, lava flows rarely cause loss of human life (Biganman et al. 2013), since their low circulation speed usually allows a timely reaction and risk prevention and mitigation measures to be put in place (Carracedo et al. 2001), such as the evacuation of people, transportable assets and belongings. That said, the effectiveness of these measures in the context of planning the response during a volcanic emergency can be more effective with improved knowledge and the ability to not only predict the behavior of the lava flow, but, above all, its likely path.
Knowledge of how lava flows move has been extensively studied by the interpretation of lava flow morphology (Hagiwara 1941;Manakami 1951;Shaw et al. 1968, Pinkerton andSparks 1978). These studies confirm the similarity of their behavior with a Bingham liquid where the flow characteristics are determined by a yield stress (Hulme 1974). The yield stress depends in turn on the viscosity of the liquid, and this depends on factors such as temperature. The physical properties, the volume of lava emissions and local conditions such as topography, are also parameters that have been studied to gain a better understanding of lava flow behavior.
In recent decades, not only have advances been made in the knowledge of lava flow characteristics, but computer science and GIS have also been improved. Thanks to this progress much of the research effort in this field has been directed towards the development of dynamic simulation models. These models try to imitate the behaviour of this phenomenon with the aim of improving both our knowledge and our ability to predict it. In general, simulation experiences are based on numerical models applied to geographic space, with a predominant use of cellular automata, and where the cell states are usually binary (lava, non-lava), e.g. Felpeto et al. (2007) although examples of prototypes that add states of solidified and fluid lava, in three-dimensional systems, were found in the literature search (Harris and Rowland 2001;Hidaka et al. 2005). In addition, there are deterministic and probabilistic approaches which specify a set of rules to drive the cell state change, e.g. nonlava to lava. The first approach provides results on discrete values and locations based on the interaction of several lava properties, such as effusion rate, temperature, viscosity, rheological properties (density, viscosity, yield stress), chemical composition, crystallinity, etc. and environmental properties, where the topography is always included as a factor (Harris and Rowland 2001;Hidaka et al. 2005;Ishihara et al. 1990;Young and Wadge 1990;Bilotta et al. 2012;Mossoux et al. 2016); the second approach estimates the most probable path and the extension of the flows by basically using environmental characteristics (relief: aspect and slope) where a set of parameters need to be calibrated, such as lava thickness or constant values to simulate the yield stress (Felpeto et al. 2001;Favalli et al. 2005;Damiani et al. 2006). These parameters can be empirically calibrated using a statistical method or by applying a machine learning algorithm.
Generally speaking, the deterministic model approaches provide accurate simulation results that can even predict the velocity, the shape and other dimensional and structural characteristics of the lava flow other than the likely path of the lava flows. Nonetheless, they usually require a set of data that is not always available in real time during the eruptive process. Moreover, during such situations the operating response teams do not need such accurate information about the lava flow evolution and its properties, but having a map showing the geographic spaces that will be invaded by the lava could be critical. From this perspective, probabilistic methods such as the LAVA FLOWS model [8] or the DOWN-FLOW model by Favalli et al. (2005) seem to be more effective in terms of the speed required during the emergency. Furthermore, these algorithms are able to run and produce predictions only with a suitable Digital Elevation Model.
In the case of the disaster caused by the strombolian-Hawaiian volcanic eruption in Cumbre Vieja on the island La Palma in the last four months of 2021, industrial buildings, infrastructures, crops, villages, and population were not only threatened by the lava flows but also affected by the different volcanic hazards. In this respect, the lava destroyed many of these assets, dramatically burying all that it found in its path. The civil protection deployment coordinated by PEVOLCA (Volcanic Emergency Plan for the Canary Islands) needed the best available information on the evolution of the volcano activity on a daily basis to manage the situation. The successive vents of the eruptive fissure were located about 900 m above an area of dispersed residential land use, small villages and a set of parallel and vital roads connecting some of the main economic activity nodes of the island. The geomorphology of the land corresponds to that of a young geological ramp of less than 20,000 years (Carracedo et al. 2001) where the topography does not always clearly show the slope and the drainage channels are still not well built (Romero 1990). In this context, PEVOLCA and its human resources (formed by members of the Instituto Geográfico Nacional (IGN), Instituto Volcanológico de Canarias (INVOLCAN), Instituto Geológico y Minero de España (IGME), Instituto Español de Oceanografía (IEO), Universidad de La Laguna (ULL) and Universidad de Las Palmas de Gran Canaria (ULPGC)) requested continuous information about the evolution of the lava flows and their likely occupation. Therefore, a probabilistic lava flow predictive model was used to provide the necessary information for risk reduction and management. In order to do this, the authors participated collaboratively, in collaboration with the Island Council technicians, to provide data from the quasi real time flow simulations and to contribute to the emergency management to reduce the risk of disaster. However, the authors were not members of the PEVOLCA scientific team.
One of the main goals of the present work was to accurately estimate the parameters of the LAVA FLOWS simulation model (Felpeto et al. 2001) used in the Q-LavHA plugin for QGIS software (Mossoux et al. 2016) to simulate the lava flows of past volcanic events in La Palma, and the application of the model in the recent eruption which occurred in Cumbre Vieja, La Palma in 2021. The use of the algorithm during the 2021 volcanic emergency, the simulation process and the technical steps followed are described below to show the effectiveness and usefulness of the model in terms of accuracy, but to also demonstrate its value as an accessible and easy-to-use tool for response management in eruptive events with characteristics similar to those in the eruptive process on the island of La Palma in 2021.

Volcanic disaster context
The island of La Palma is one of the seven islands belonging to the Canary archipelago and is considered the most volcanically active island, given that seven of the fifteen historical eruptive events that have taken place in the Canary Islands in the last 500 years, since the first records began, have occurred in La Palma (Troll and Carracedo 2016). La Palma is an island with a triangular base, which can be differentiated into the following two large topographical units: the Taburiente Volcano, in the north, and Cumbre Vieja, in the south. The latter is where volcanic activity has been concentrated for 150,000 years (De la Nuez et al. 2018). Cumbre Vieja is a volcanic ridge characterized by monogenic eruptions which occupies an area of 220 km 2 and has a length of 20 km arranged in a north-south direction (Fig. 1a). Although emission centers can be found distributed throughout its surface, with poorly defined fractures, there was a reorganization of volcanism twenty thousand years ago and these fractures began to line up in the center of the ridge, in the current NS axis (Day et al. 1999). All the historical eruptions that have taken place on the island are classified as being fissural in nature, in other words, they have the appearance of several emission centers distributed along eruptive alignments that function simultaneously or alternatively. In addition, all the flows from these eruptions, with the exception of those from Montaña Quemada (ca. 1430-1447), have reached the sea and increased the surface area of the island due to the formation of lava deltas. It should be noted that the geologically young age of the area means that there are few erosive forms. Even so, there are indications of an incipient drainage network with simple channels, without defined headwaters, or even without the channels being connected to the sea (Romero 1990).
In this context, on September 19, 2021, a new eruption began in the northwest sector of Cumbre Vieja and lasted for eighty-five days, becoming the longest recordedhistorical eruption on the island and breaking the trend that these were of increasingly shorter duration (Troll and Carracedo 2016). In line with the previous historical eruptions in the Canary Islands, the volcanic event was characterized being of a strombolian fissural eruptive typology with phreatomagmatic pulses, reaching Volcanic Explosivity Index (VEI) of 3. In addition, the fissural nature of the eruption led to the opening of eleven vents over the course of the eruptive process, some of which expelled lava simultaneously, a fact that caused the overlapping of lava flows at times. As a result, after the end of the eruption, six craters were aligned in a NW-SE direction along a 557 m fissure. On the other hand, according to the data provided by the PEVOLCA scientific committee, the total amount of material emitted was more than 200 Mm 3 , and the estimated area occupied by the lava flows was 1200 ha. Similarly, the lava flows, which are characterized by the majority presence of aa-type lavas, have a mean average thickness of 12 m that, on occasions, reached maximum heights of up to 70 m. Besides, with regard to the subaerial route of the lava, its longitudinal development is 6.5 km and its cross section is 3.3 km at the most equidistant points. Furthermore, some of these lava flows entered the sea, giving rise to an increase in the island's surface area of 48 ha by creating two platforms in the form of fans, known as lava deltas. These lava structures, in turn, are supported by the materials from the emitting centers whose advance below the water is up to 1.1 km, occupying an area of more than 21 ha (IGME 2022).

Data source
The present work uses sources of information on the geographic data necessary for the simulation of lava flows. The data used are mainly divided into three groups: elevation models, location of emission points and evolution of the flow perimeter. Two different Digital Elevation Models (DEM) were used over the course of the eruptive process, both with a resolution of 5 m. The first DEM used was prior to the start of the eruption and is available in the open data of Instituto Geográfico Nacional (IGN), the institution in charge of volcanic surveillance in Spain. A new DEM, including the changes in the previous topography caused by the eruptive process, was built during the course of the eruption. In addition, the perimeters of the lava flows shared on the open data website 1 3 of the Cabildo de La Palma (https:// www. opend atala palma. es/) (Cabildo de La Palma 2022) have been used to perform the validation of the simulations.
The location of the emission centers was initially obtained from the observations of the local emergency teams deployed on the ground. Subsequently, the activation of the remote imaging services provided by the Copernicus program allowed a more precise location of the emission centers.

Methods
Given the objectives of the present work, the methodology described below is not only aimed at describing the model used and the validity of the results, but also presents the experience of lava flow simulations before and during the volcanic eruption in Cumbre Vieja in La Palma. In this regard, (1) the choice of the model used, its operation and the justification of its previous calibration results are described here (Martín-Raya 2020); (2) data were exchanged between the technicians belonging to the emergency groups and the research group to perform the simulations; and (3) the method and results of the validation of the simulations performed during the eruption are explained.

Model selection, operation, previous calibration and validation
Given the urgent need for information from PEVOLCA and considering the difficulties concerning the availability of all the parameters related to the characteristics of the lava (rheology, viscosity, temperature…) and to the eruption itself (such as the rate of emission or the flow rate), the use of deterministic models was ruled out and a probabilistic algorithm with low input data requirements was chosen instead. In this case, the LAVA FLOWS probabilistic model (Felpeto et al. 2001(Felpeto et al. , 2007 was selected, which requires few minimum parameters and, therefore, is of great value and utility in emergency situations. This model simulates the probable trajectories that the lava flow fronts will take assuming that the topography plays the main role in the advance of the flow. In addition, this model is incorporated into the Q-LavHA software (Mossoux et al. 2016), a QGIS plugin that is both easy to use and readily available. Felpeto's model (Felpeto et al. 2001) expresses, from a digital elevation model (DEM), the probability that a cell is flooded by lava. An algorithm based on the Monte Carlo method is used for the calculation of probability. Starting from the emission center, the model assumes that the flow propagates from the central cell towards one of the eight cells surrounding it, calculating the flow probability (Fig. 2), based on the altitude differences, by applying the following formula: where P i is the flow probability, Δh i is the difference in height between the central cell and each of the eight surrounding cells.
It should be remembered that if the difference ( Δh i ) is negative, the probability will be 0, because the propagation of lava uphill is impossible. However, the model has corrective factors that allow the lava to flow over small obstacles or depressions, so the thickness ( h c ) of the lava is always added to the value of each cell before calculating Δh i . Furthermore, since lava can fill in pre-existing depressions and continue advancing, a larger corrective factor ( h p ) is introduced. This is applied if the cell is completely surrounded by others of a greater height, and even considering the height of the lava flow, the lava flow could not flow over them. Therefore, as can be seen in Fig. 3, the simulation progresses smoothly until it encounters a depression that it is unable to move over. Therefore, the h p correction is added and when it exceeds the next pixel elevation, the simulation continues.
After determining which cells could be affected by the flow by calculating the probability for each one, a random cumulative probabilistic value ( S i ) is calculated by applying the Monte Carlo method to select the cell which the flow will propagate across. Several paths are calculated where the flow cannot return backwards and the highest probabilities coincide with the areas with the greatest slopes. After each path, a probability of being affected is laid on each pixel ( P j ). This probability is added with the rest to define the final pixel probability. The described model allows the configuration of the desired number of iterations when performing a simulation, although it should be considered that the recommended minimum to obtain good results is 1500 (Mossoux et al. 2020). After this process, and following Monte Carlo method, 1 3 a random number between zero and one is drawn by the algorithm and compared with S i . If the random number falls within the interval [ S i − 1 , S i ], the pixel is selected as the next pixel on the lava flow path (Mossoux et al. 2016). The cell will be flooded by the lava flow if: P j : probability of being affected in each iteration using Monte Carlo.n rnd = random number between 0 and 1, calculated by the algorithm each time.
If the result is negative, this means that the lava has entered a depression that even when applying the corrective factors, it is not able to flow over it and, therefore, the model stops.
Knowing the final location of the lava flow and determining at what point it will stop is not a straightforward process, since it depends on complex factors, which are difficult to obtain and vary among eruptions (Daminani et al. 2006). However, the Euclidean maximum length (L max ) option incorporated in the model itself was used to prepare a quick and useful forecast for the emergency situation. This value can be easily estimated by studying, for example, the lengths of historical lava flows. When the lava flow reaches the set distance, the model ends. For this reason, and given the characteristics of this specific case, the simulations performed here are parameterized to force the lava flows to reach the sea, given that the main demand for information involves the trajectory of the fronts and their foreseeable occupation.
Q-LavHA has an algorithm that, by means of a fitness index, allows the simulation of the lava flows to be validated by comparing them with a real lava flow, similar in origin, composition and morphology. Therefore, a previous step to performing the simulations was the calibration of the parameters used to achieve a good goodness of fit through an empirical trial-error method (Table 1). Here, the lava flow from the eruption of the San Juan Volcano in La Palma (1949) was used to perform the calibration (Fig. 4). The optimal way to carry out this calibration would be to use the pre-eruptive topography as is usually done (e.g. Rodríguez-González et al. 2019;Becerril et al. 2021). However, the first available topographies of the island of La Palma date from the late 1960s and there are no previous elevation models. As such, there is also no information on the thickness of the flows, which is needed to accurately eliminate the lava body. For both of these reasons, the DEM for the year 2021 was used here.
In addition, data on the historical eruptions in the Canary Islands were considered, such as the thickness of the lava flow (Carracedo et al. 2004) and the maximum length of the lava flows in eruptions in La Palma. Regarding the latter, the greatest length reached by a lava flow in the volcanic eruptive events in La Palma is ten kilometres in the case of the Montaña Quemada eruption (ca. 1430 − 1447). However, although most of the flows have maximum visible lengths notably shorter than those of Montaña Quemada, with paths that do not exceed eight kilometers, this distance (10 km) was considered appropriate to ensure that the model simulated lava flows which are able to reach the coast, as happened in six of the seven historical eruptions on La Palma (Romero 1990). The fitness index compares the overlapping (true positive) between the real flow and the simulation. After conducting numerous tests, the highest fitness index value obtained was 0.3 (Fig. 4), that is, 30% of the simulation adjust to the real flow. It should be said that this value would improve if the pre-eruptive topography were used. Although Mossoux et al. (2016) suggest 0.5 as the optimal value, 0.3 is considered sufficient due to the objective of the calibration. Furthermore, false positives and negatives are also an important consideration to validate the calibration of the parameters, so these values must be taken into account to correctly interpret the results. For effective use in risk management, overestimating (false positive) implies applying an acceptable precaution in the context of hazard assessment, while underestimating (false negative) the affected area would be problematic (Mossoux et al. 2016). In this regard, although 70% of the surface is overestimated (0.7 false positive), only 5% is underestimated (0.05 false negative), which means it is useful in an emergency.

Development of simulations during the volcanic emergency
The objective of the lava flow simulations in the present work is to provide information for the emergency management teams about the possible trajectory of the flows, especially during their progress on the ground. As explained above, the geographic location of the volcanic event corresponds to an initial location of an eruptive fissure at high levels on the island (900 m) in a location with a small population and little human activity, but from which lava flows run towards the coast (elevation 0). In this context, emergency managers and citizens require information quickly about the possible trajectory of the emerging and advancing lava flows that are destroying buildings, residential neighbourhoods, crops, and other infrastructure at lower levels which ends up buried under the path of the lava as it moves towards the sea. However, the direction of the lava flows cannot be easily predicted because the channel networks are not clearly defined due to the geologically young age of the area and its little erosion.
As the volcanic episode unfolded, the data and the flow of information generally improved, as a result of the work of the PEVOLCA scientific committee teams and the observation of local emergency operations, as well as national detachments specialized in the response to disaster situations, such as the Military Emergency Unit (UME, in Spanish), in addition to remote observation mainly provided by the European emergency and early warning service Copernicus. There were constant changes in the behaviour, evolution, extension, and displacement of the flows, as well as the emission points and the morphology of the vents. Therefore, the data and information generated were also highly dynamic. Furthermore, and along the same lines, the demand for predictive information on the progress of the lava flows became a critical factor in monitoring the entire volcanic process.
In short, the entire eruptive episode and the changes in the evolution of the flows, the migration of the eruptive vents and the changes in the topography themselves, varied as much as the availability of data and information does. In addition, the simulations carried out can be classified according to the evolutionary circumstances of this specific eruptive event and the availability of data during it, and as such the following simulations are described below: Simulation 1 This was performed at the beginning of the volcanic episode between September 19th and 31st 2021, when the cone had not yet been defined and, although there were three main vents, up to nine emission points were counted along the fissure. The lava flows were aa type, highly viscous, with mean average speeds of 200 m per hour, decreasing to less than 1 m/h in some stages, especially in the furthest points from the eruptive vents (IGME 2022), and notable thicknesses were observed in the direct observation of lava flows that, in some of the more depressed areas, exceeded fifteen meters. The first simulation was performed on September 20th, using both the imprecise location of the aforementioned three main vents and the parameters of the previous calibration (Fig. 3) as well as the DEM available in the abovementioned open data of the IGN.
The algorithm in Q-LavHA allows the configuration of the model to simulate lava flow emission either from all three vents at the same time, from a series of points evenly arranged along a fissure or from each vent individually. The last option was chosen, subsequently adding the probabilities of being flooded in each simulation. Thus, the areas that were in danger of being occupied by the lava flow were identified and the probabilities of such an event were obtained. To better understand the existing danger in each area, the resulting probabilities were defined using the quantil method as this has been used in other cases with other hazards (Tehrany et al. 2014;Jaafari et al. 2014;Rahmati et al. 2016) (Fig. 6), and with extremely low values being discarded (Table 2).
Simulation 2 This simulation corresponds to the intermediate stage of the event between October 1st and 20th. During this period, the previous lava flows, as well as the fall of pyroclasts and ash, modified the topography of the affected areas. At the same time, the outer layers of the flows cooled down and the formation of volcanic tubes and lava channels between levees occurred. As a result, the lava flows, which were generally more fluid than in the initial stage, mainly circulated on top of the lava flows emplaced in the first stage of the eruption. However, the appearance of a new emission point, to the north, needed new simulations to predict the paths of the lava flows that threatened to cover new land. On the other hand, the main front of the initial lava flows had already reached the sea, but new lava flow fronts, generated by processes of new lava flowing over existing lava, threatened new areas which had not yet been affected.
The location of the new eruptive vent located further north and the perimeter occupied by the lava flows were required to perform this simulation. In addition, it was necessary to obtain information on the new topography generated by the materials ejected by the volcano, mainly the lava flows. Given the rush and lack of new terrain elevation data, a decision was taken to perform a simulated mean average elevation of the terrain occupied by the previous flows using a calculation based on observation and data taken in the field work. In this respect, a mean average extrusion of seven meters was made to all the cells occupied by the polygon representing the extension of the flows at the time. This extrusion contributes, above all, to modelling the effect of the levees and fronts of the previous flows in a more realistic way. These morphologies form important obstacles to be overcome by the active flows running adjacent to or with transversal paths to the previously partially or totally cooled flows.
Simulation 3 This simulation concerns the last stage which was specifically studied for this eruption and, above all, in the context of the present work, covering the period from October 20th to December 14th, 2021. As observed throughout the first month of eruption, the flow thickness was greater than the simulated thickness. As a result, it was decided to change the simulation parameters to be consistent with the real situation, by adjusting the values of H c and H p to 7 and 4 m, respectively. The lava flows generally continued to move, for the most part, either on top of the previous lava flows, or through volcanic tubes. Several lava channels were seen running to the sea without occupying new land. However, some fronts, even with unoccupied space ahead, advanced slowly. Lateral thickening of lava flows was also produced by new lava flowing over existing lava. These overflows sometimes formed lava streams that flowed alongside the larger lava flows, forming small fronts similar to those directly produced by an emission point. This situation occurred, at this stage, especially in the south since the flows from the north flowed mostly over the pre-existing flows. This is where a new vent appeared in the altitudinal strip of the initial fissure and it was therefore necessary to perform a new simulation. This vent to the south, which initially only expelled gases, ash and some pyroclasts, began to extrude lava around October 25th, affecting new areas in the southernmost part, burying part of a cemetery at around 550 m above sea level. Below this altitudinal strip, the advance of the lava flow to the south was mainly produced by the abovementioned overflows (Fig. 5). Although consideration was given to simulate these overflows which were taking place at lower levels, it was decided to use the modelled flow for the southern vent as a reference. To a large degree, this simulation also covered the terrain that could be occupied by lava from the aforementioned frontal and lateral overflows at lower levels.

Validation method
Since the main objective of these simulations is focused on predicting the trajectory and occupation of the flows, the index proposed by Cohen (1960) called the Kappa Index was used to perform the validation of the simulations. This is a measure that expresses the coherence existing between objects observed in a two-dimensional plane, in this case, between the lava flow simulations and the actual lava flows of the eruption. Its computation focuses on the proportion of agreement between the two compared objects, taking into account the overestimates and underestimates of p c , in this case between two cell states (lava or non-lava), and is calculated as follows (Walker 2010).
where p c corresponds to the expected corrected proportion and p 0 is the observed corrected proportion. These parameters are calculated as follows: where TP is true positive from the simulated lava flow, TN: true negative, FP: false positive, FN: false negative and N is the number of pixels in the map (Fig. 6) (Sahani and Gosh 2021). The results of the calculation are between 0 and 1, where 1 corresponds to perfect consistency.
The Kappa index analyses and compares the cells of each of the layers, establishing where the agreements and discrepancies are found (Fig. 5). To correctly apply the Kappa index in a geographic information system, it is necessary to standardize the information to be used. Therefore, the information must be in raster format, have the same extension and be classified in values 1 and 0 (where 1 is lava and 0 is non-lava). It should be noted that, in this case, to make the comparison between the real flows and the simulations, the decision was taken to eliminate the new lava deltas formed offshore, since the DEMs used in the model do not contain bathymetry data. In the first moments of the eruption, the simulations occupied spaces that the real lava flow did not, and which were later filled with the advance of the lava. Therefore, in order to determine whether the simulation was accurate or not by using the Kappa index, it is necessary to wait for the end of the eruption when the lava has already occupied the whole territory. However, it was decided to subject the simulations to a successive comparison of the simulation results with the area occupied by the lava flows at successive times of between five and ten days, depending on the observed rate of volcanic activity. This makes it possible to observe how the Kappa index approaches 1 (greater degree of similarity) as the perimeter of the land occupied by new flows increases or, in other words, to check the validity of the simulations over time.

Results and discussion
The results presented below are mainly aimed at describing the simulation process and the technical steps followed, as well as the capacity of the model to predict the path of the lava flows using a validation procedure based on the Kappa index.
As described in the methodological section above, based on a series of significant changes in the evolution of the lava flows within the temporal evolution of the eruptive phenomenon, three simulations were required at different times during the event. The cartographic results are shown in Fig. 7 (see dynamic version for Fig. 6).
The simulations were performed assuming that the volume of lava emitted was high enough for the lava front to reach the coast, and as such the simulations are considered valid when the Kappa index tends to reach the value 1 (maximum similarity) (Fig. 7).
During the first half of the eruption, where the lava repeatedly occupied new ground, the index grew rapidly, from 0.212 to 0.731, on October 30th (Fig. 8, day 41). After this, a slower increase of the index was observed until it reached the maximum of 0.793. This slowdown coincides with the stagnation of the affected surface, whose evolution is similar. During the first forty-one days of the process, 965 ha were occupied, almost 80% of the Fig. 6 Schema of Kappa index. It compares simulated and real lava flow occupation and estimates the similarity between each of them total, as shown in Fig. 9, while in the second half of the eruption process only 255 more ha were invaded (the remaining 20%). This is because the lava began to flow more fluidly over already covered land or through volcanic tubes.
The noteworthy falls in the goodness of fit produced by the Kappa index mainly occurred in the intermediate stage of the eruptive process and there was also a less pronounced one in the final stage (Fig. 9). The falls in the Kappa index were produced by the addition of a new simulated lava flow and, in both cases, this was followed by an increase in the Kappa index because the lava flows covered the predicted area, which reinforces and strengthens the accuracy of the model. On the other hand, the slight decrease that occurred at the end of the process is related to the overflows that occurred around the cone as a result of an intensification of the eruption due to the final degassing process. This becomes clear by observing the evolution of occupied hectares, where a more pronounced increase in the trend was observed in the final two days of the eruptive process, occupying 4% of the total area affected, that is, 41 ha. This explains why, although an index value of 0.8 was reached on December 12th, the final result is slightly lower (0.793).
It should be noted that two discordant zones stand out at the end of the eruptive process: one spatially overestimating and the other underestimating the area. The first corresponds to the southern strip (simulation 3 in Fig. 7), where the short duration of the eruption left the areas of moderate and low susceptibility unoccupied, as mentioned above. On the other hand, around the cone, a marked discordance can be seen, due to the continuous changes in its morphology, related to the process of construction and reconfiguration. However, according to the classifications for interpreting values derived from the Kappa index prepared by Landis and Koch (1977), Fleiss (1971) or Cicchetti (1984), the results obtained are relevant and adjust with a high level of accuracy to the reality observed during the entire eruptive event. In addition, the success and usefulness of the simulations require knowing the precise location of the emission points as soon as possible (Capello et al. 2016) to favour effective response processes. However, providing accuracy rapidly is a complex matter as this emergency has shown. The maximum accuracy of the location of the vents can be dispensed with to a certain extent, especially at the beginning of the episode bearing in mind the complexity of fissural eruptions in terms of the location of emission points.
Ultimately, this work not only validates and complements, but also shows the high applicability of the previous studies of Felpeto et al. (2001Felpeto et al. ( , 2007. In the sense that lava flows are complex systems whose behavior and evolution are driven by the interaction of many factors such as temperature, degree of crystallization, rheology and the environment (Bilotta et al. 2012). Many of the models that simulate the spatial and dynamic behaviour of a natural or human phenomenon using cellular automata are based on the complexity theory (Díaz Pacheco 2015). This means that, among other things, these models are able to find patterns in apparent chaos and imitate the interactions between factors and elements (temperature, topography, chemical composition…), reproducing, in this case, the behavior of the lava flows using a set of simple rules (Tobler 1979;Couclelis 1985;White and Engelen 1993). When information from past eruptions exist, it is possible to compile data about those parameters and use them to simulate lava flows (Harris et al. 2011;Zuccarello et al. Fig. 7 Results of lava flow simulations run during the volcanic eruption. The first two runs were performed using the DEM available before the eruption and the last one with the new updated DEM. Source: IGN (DEM 1); Island Council (DEM 2,3; lava perimeters; fabric) 1 3 2022) or alternatively real-time measurements of time-averaged discharge rates derived from satellite thermal infrared images can be used (e.g. Vicari et al. 2011;Kaneko et al. 2022) or even data from other eruptions in similar volcanic contexts can be used. These deterministic models usually produce a more accurate prediction of the likely paths and extension of lava flows, especially when the vents are located far from the coastline and the final length of the flow is uncertain. However, in a complex and sudden emergency situation, as was the case in La Palma, it is necessary to know as quickly as possible which areas would be affected. Therefore, a probabilistic model, despite being less accurate, is shown here to be the better option where previous data did not exist and real-time data are poor, particularly at the beginning of the eruption. One of the main advantages of the method used in the present work is its ability to simulate a complex natural event and/or a non-linear phenomenon using a set of simple rules integrated into one simple equation.
Many of the scientific works on lava flow simulation have traditionally been focused on the characterization of the physical phenomenon (Sobradelo et al. 2011;Barker et al. 2015;Branca and Del Carlo 2005;Cronin et al. 2003;Umakoshi et al. 2001) but there are few works focused on the value of lava flow simulation as a risk management tool (e.g. Favalli et al. 2009Favalli et al. , 2012Becerril et al. 2014). Similarly, the results obtained here show the great value of the method designed here as an early warning system for monitoring and forecasting a threat such as lava flows. In addition, elements such as the predictive cartography developed during the eruptive process in La Palma as a result of the continuous simulations helped to specify the decision-making associated with a danger such as lava flows. The latter is in line with that expressed in international reports such as the Global Assessment Report (2019), which emphasizes the importance of taking actions to improve the monitoring, analysis and forecasting of threats and their possible consequences. This, among other things, has made it possible to enhance emergency management processes and guarantee correct decision-making in highly relevant aspects such as preventive evacuation and multiple logistics tasks related to emergencies (Marrero et al. 2019).

Conclusions
The first steps of the present work were based on the calibration of the necessary parameters to use the LAVA FLOWS algorithm (Felpeto et al. 2001) during the development of the Tajogaite volcanic episode in 2021 on the island of La Palma. In order to do this, data from historical eruptions that occurred on the island since 1400 were used, especially those concerning the San Juan eruption in 1949, since part of its lava flow is close to the 2021 eruption, showing similar lava flow behavior. This previous calibration allowed our research group to use the algorithms to periodically supply early warning information to the operative emergency team deployed in La Palma. The authors of the present work did not have any role in the PEVOLCA (the regional emergency plan for volcanic risk) Fig. 9 Temporal evolution of the Kappa index (principal axis) and hectares of land covered by lava (secondary axis) during the eruption. Source: Island Council (lava surface) scientific committee, but they voluntarily supplied this kind of information, which was requested by the island government emergency operation team. In this regard, this work attempts to show the results of a collaborative experience during a volcanic emergency, between academics and the local administration. Everything functioned as a coordinated and complementary support to the official work of PEVOLCA.
The algorithms used are incorporated in the Q-LavHA plugin (Mossoux et al. 2016), which is easy to install and use in the free software QGIS, which makes it an accessible tool, quick to use and readily applicable during an eruptive process requiring rapid and precise predictive information about the evolution of the flows, especially in terms of their possible trajectories and their occupation. In fact, this work also evaluates the capacity of such algorithms implemented in applications of easy availability to be used in response situations to a disaster of volcanic origin.
The results show that the model used here is of great value in meeting the demand for information required by emergency operations. Its versatility, efficiency and low input data requirement provide good simulation results during the risk management tasks to mitigate the dangers inherent during a volcanic event and it does this in practically real time. At least, and the present work demonstrates this, for the type of volcanic events with characteristics similar to the one that occurred in La Palma in 2021.
The following three points should be mentioned, which could easily be improved for future experiences. The first is that the accuracy of the calibration fitness index could be improved by subtracting the topography generated by the lava flows of historical eruptions used for calibration; secondly, it is crucial to have a more accurate location of the active vents to improve the lava flow prediction; and thirdly, it is important to mention that, after obtaining a low fitness index and a reasonable location of the vents, the results of the simulations during the emergency showed a high degree of accuracy. This latter point may be due to the regular topography of the area where there are ravines or canyons which are not particularly abrupt.
Reducing the impact of a disaster requires anticipation and the present work also highlights the importance of having systems in place capable of generating accurate information during the occurrence of a volcanic event, especially at the beginning of the eruption. Information about the location and the evolution of the eruptive vents is critical to simulate lava flows and to use this type of probabilistic model based on cellular automata. An accurate DEM (5 × 5 m is good enough), and a system which is able to update this model in real time are also indispensable elements for lava flow simulation. In addition, as this experience also shows, in order to successfully channel this type of task, aimed at improving early warning systems, it is important to maintain constant feedback between the observation teams deployed in the field and the agents and scientists who perform the simulations. declare that they have no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.