Modelling and simulation of a lava flow affecting a shore platform: a case study of Montaña de Aguarijo eruption, El Hierro (Canary Islands, Spain)

ABSTRACT Recent subaerial volcanism at El Hierro Island (Canary Islands, Spain) consists of monogenetic volcanic fields. This volcanism generated cinder cones, tephra air-fall deposits, and lava flows. The lava flows reach several kilometres in length extending through shore platforms and, sometimes, penetrating under the sea level. The volcanic landforms of El Hierro convert it into a natural laboratory for topographic and morphometric modelling and lava flow simulations. We perform the modelling and simulation of the Montaña de Aguarijo eruption, a cinder cone at the NE rift. The associated lava flow channelled through a V-shaped ravine until reaching a cliff, where formed cascades. The flow spread at the cliff base over a platform before reaching the sea modifying the coastline. Different maps were designed to show the results, including the geomorphologic reconstruction of the area affected by this eruption and the lava flow simulations obtained with the Q-LavHA plugin.


Introduction
Modelling and simulating volcanic lava flows help evaluate volcanic hazards associated with monogenetic basaltic eruptions and support safe strategies. The resulting maps are critical when applied to volcanic zones where lava flows represent the main hazard for people and properties (e.g. Nieto-Torres et al., 2021;Rodriguez-Gonzalez et al., 2021). While some maps are based on qualitative analysis of historical eruptions, modern assessment combines the geological history with probabilistic or deterministic computational models (e.g. Cappello et al., 2016;Cordonnier et al., 2015;Damiani et al., 2006;de' Michieli Vitturi & Tarquini, 2018;Dietterich et al., 2017;Favalli et al., 2009;Richter et al., 2016;Sieron et al., 2019;Tarquini et al., 2019). The Digital Elevation Model (DEM) of pre-flow topography and the vent location are essential. In addition, physical and rheological parameters of magma are required for simulation using a deterministic approach. In this work, preand post-eruptive DEMs are the product of an accurate 3D palaeogeomorphological reconstruction of landforms before and after the studied eruptive event, based on fieldwork. In both cases, the tridimensional reconstruction of palaeosurface follows the methodology of Rodriguez-Gonzalez et al. (2010), aiming to be closer to reality than previous formulae-derived methods.
Fieldwork is the first stage to obtain the necessary information for reconstructing the topography and morphology of past eruptions from which the preand post-eruptive DEMs are derived. It consists of detailed geological mapping of the main volcanic units (cones, tephra air-fall deposits, and lava flows), together with the acquisition of information of geomorphologic (e.g. lava flow thickness), structural, stratigraphic, and, if necessary, bathymetric observations ( Figure 1). This work aims to provide an effective workflow for the pre-and post-eruption topographic and morphometric reconstruction of monogenetic volcanic eruptions and assess related lava flow simulations ( Figure  1). As a case study, we test its effectiveness in Montaña de Aguarijo, one of the recent and well-preserved monogenetic volcanoes of El Hierro Island (Canary Islands, Spain). The eruption of Montaña de Aguarijo has the peculiarity of exhibiting a lava flow forming cascades on a coastal cliff and spreading at the cliff base over an old coastal platform reaching into the sea. main islands and several islets, are located off the north-western African coast, between 29°25' and 27°37' N and 18°10' and 13°20' W ( Figure  2a). This archipelago is preceded by a group of seamounts located further NE and together form the Canary Volcanic Province (CVP), an ∼800-kmlong and ∼400-km-wide volcanic belt. The CVP developed above the Jurassic ocean lithosphere (cold, thick, rigid, and old crust) lying close to the African passive continental margin and on the slow-moving (∼2 cm/y) African plate (Silver et al., 1998). Volcanism in the CVP decreases in age from the NE to the SW, from ca. 68 Ma in the Lars/Essaouira seamount (Geldmacher et al., 2005) to ca. 1.12 Ma in the El Hierro island (Guillou et al., 1996) (see Figure 2a). This age progression is due to the African plate's movement over a mantle plume (e.g. Carracedo, 1999;Geldmacher et al., 2005;Hoernle & Schmincke, 1993;Zaczek et al., 2015).
The Canarian archipelago displays, similarly to other intraplate volcanic islands, three main hotspot volcanic stages of evolution: (1) juvenile (shield) stage, (2) volcanic quiescence stage and (3) rejuvenated stage. The westernmost islands (El Hierro and La Palma) are in the juvenile stage, which is   (Geldmacher et al., 2001;Guillou et al., 2004;Ryan et al., 2009). (b) El Hierro Island ortophotomap and location of Montaña de Aguarijo. characterised by high eruption rates and fast volcanic growth, whilst the easternmost ones (Fuerteventura and Lanzarote) are in the rejuvenated stage where the erosive processes play a significant role (Carracedo et al., 1998(Carracedo et al., , 2002Walker, 1990).

El Hierro
El Hierro, the youngest and westernmost island of the Canarian archipelago, is the emergent area of a 280 km 2 volcanic shield that rises from ∼4000 m deep seafloor to 1502 m above sea level at the centre of the island (Pico Malpaso) (Figure 2b). The most characteristic feature of El Hierro is its tetrahedron shape, with edges formed by three convergent ridges of volcanic cones at 120°and separated by wide embayments.
Montaña de Aguarijo is a cinder cone located at the north-eastern rift system within 3 km to the capital Valverde and about 5 km to the small airport (El Hierro Airport). Together with Montaña del Tesoro and Montaña Chamuscada represent the principal platform forming eruptions at the NE rift and are thought to be younger than 18 ka BP due to older lavas are arranged to form cliffs (Carracedo et al., 2001). The most particular feature of Montaña de Aguarijo is its 2300 m length lava flow that cascades over the edge of a palaeocliff and spreads over a narrow old shore platform before reaching and plunging into the sea (Figure 3).

Mapping and palaeogeomorphological reconstruction
Volcanic morphometric modelling included intensive fieldwork to carry out a meticulous geomorphological and topographic reconstruction of the Montaña de Aguarijo eruption, to obtain the pre-and post-eruption DEMs. Fieldwork identified and analysed the volcanic landforms, i.e. cone, lava flow and levées, besides the underlying palaeorelief (Rodriguez- Gonzalez et al., 2010Gonzalez et al., , 2012. The geological mapping of these units was performed through a digitising tablet (iPad Pro 11" Wi-Fi + Cellular model + GPS / GNSS receiver) to trace over and vectorise the contact lines and geolocation points (photos, sketches, samples, etc.). For more accurate results, we try the heads-up digitising method with pencil, by which the operator is continuously geolocated, using the FieldMove application. We use high-resolution raster maps as MBtiles files as a reference (background) layer for the digital base map data (e.g. orthophoto, LIDAR, topography) in the FieldMove application. Since the new vector elements are displayed over the reference map image, we can visually check the accuracy as we go along. The vector data structure allows representing the geological data in a single vector object providing the most efficient means to digitise and store the data for a large, complex map. Later, for more flexibility in the map design, we may separate the data into different thematic layers (such as unit contacts, outcropping structures, or individualised eruptions). Finally, the data acquired is exported to a KMZ file and introduced in a GIS environment, which offers valuable capabilities for estimating accurately the morphological parameters and derivates required, e.g. eruption output volume, slope and aspect of the terrain (Figure 1).
Field-derived data were combined with the available 1:5,000 topographic (Cartográfica de Canarias GRAFCAN, 2006) and 1:1,000 bathymetric (Dirección General de Costas, 2003) digital maps (5 and 1 m contour line equidistance, respectively) to manually modify the present-day contour maps to derive the preand post-eruption relief. These modified contour maps were used to obtain the corresponding DEMs in the two different stages: (1) before the eruption modifies the relief (pre-eruption DEM) and (2) just at the end of the eruption, prior to erosive processes (post-eruption DEM). The comparison between these two DEMs and the present-day landform allow us to identify the geomorphological evolution of the study area ( Figure 4). In this evolution, three coastlines are recognised (pre-and post-eruption and current coastline). In this scenario, the post-eruption and current coastlines reflect the lava flow fronts that gained land from the sea and later receded due to coastal erosion, respectively. Besides, the distribution and morphology of the bathymetric curves permits the drawing of the advance of the Montaña de Aguarijo lava lobe in a submarine regime.
Obtaining the morphometric parameters of the volcanic cone and lava flow requires different processes depending on if they have a horizontal or vertical component. Horizontal parameters (i.e. lava flow length or crater rim major and minor axes) were obtained by measuring geotools on GIS software. Vertical parameters (i.e. lava flow thickness and cone heigh) were derived from comparing the pre-and post-eruption DEMs.
The bulk erupted volume (V bulk ) was calculated through GIS methods (cut-and-fill analysis) from the elevation difference between the pre-and post-eruption DEMs. The obtained value considers the solid and void volumes of the different volcanic units and, therefore, does not represent the absolute volume of the emitted volcanic products. Correction factors assuming a porous volume fraction of 75% for the volcanic cone (Mangan & Cashman, 1996) and 25% for the lava flow (Wolfe et al., 1987) were applied to obtain the dense rock equivalent volume (V DRE ) (e.g. Rodriguez-Gonzalez et al., 2010. Comparing the post-eruption and present-day DEMs also helped understand the magnitude and distribution of the erosive processes that have affected Montaña de Aguarijo since its eruption.
The pre-eruption relief of any eruption zone constrains the lava flow velocity, propagation direction, and volcanic cone shape. The pre-eruption DEM of Montaña de Aguarijo and the obtained morphometric parameters were used as input data for the lava flow simulations.

Q-LavHA simulation models
The Q-LavHA v3.0 freeware plugin (Mossoux et al., 2016), integrated into the QGIS v3.10 software, simulates the probability of inundation by 'a'ā lava flows from one or more eruptive vents spatially distributed on a DEM. Q-LavHA integrates three different models that allow establishing a simulated lava flow path, where the first two are probabilistic and the latter deterministic: (1) Maximum Length (L max ), (2) Decreasing Probability (L normal ), and (3) FLOWGO (L flowgo ).
The L normal alternative has not been applied since it is based on a decreasing cumulative density function (equation 9 of Mossoux et al., 2016), and therefore, requires the average length and standard deviation of previous lava flows with similar behaviour. In this work, we present Montaña de Aguarijo's lava flow simulation running L max (Felpeto et al., 2001) and L flowgo (Harris & Rowland, 2001) models. In both cases, corrective factors (H c , minimum thickness and H p , mean thickness) and the length (L) of the actual lava flow are considered to determine the propagation of the simulated flow and allow it to overcome small topographic obstacles and fill depressions. The number of iterations is also specified, which define the number of flow lines computed to obtain a good fit between the simulations and the actual lava flow (Mossoux et al., 2016).
The probabilistic maximum length model (L max ) considers a maximum length (L') as the distance covered by the lava flow line, estimated by applying a corrective factor to the actual lava flow length. The simulated flow runs from a given emission point over a DEM based on propagation rules, and each iteration stops when the lava flow line reaches the specified maximum length (Mossoux et al., 2016). The deterministic FLOWGO model (L flowgo ) requires the input of multiple physical and rheological parameters that determine the evolution of the simulated lava flow over the topography. The simulation stops when at least one of the following conditions is achieved: (1) the lava velocity is zero, (2) the lava core temperature reaches the solidus, or (3) the yield strength at the base of the channel is greater than the downhill stress (Mossoux et al., 2016). The slope is a critical factor in FLOWGO since it influences the lava flow propagation and cooling speeds, two key parameters that control the final length reached by the lava (Rodriguez-Gonzalez et al., 2019).
To evaluate the accuracy of each simulated lava flow, six fitness indexes (FI; whose value ranges from 0 to 1) are calculated by Q-LavHA through the same DEM as in the simulation but assigning a value of 1 to the pixels that represent the actual lava flow and a value of 0 to the pixels outside of it. Fitness indexes determine the overlapped (FI true positive , TP), overestimated (FI false positive , FP) and underestimated (FI false negative , FN) areas between the simulated and the actual lava flow. Therefore, the closer the TP value is to 1, the better the correspondence is between both flows. The FP and FN indexes assess whether the divergence between the simulated and actual lava flow is due to an overestimation or underestimation of the area covered by the simulation, helping determine its effective use in hazard management properly. A weighted composite score (FI composite score , CS) is also obtained from the TP, FP and FN indexes, whose value reflects the null (0) or maximum (1) coincidence between the simulated and the actual lava flow. Besides, the cumulative flood probability (FI cummulative probability , CUM) is established for the pixels located inside (FI CUMin ) and outside (FI CUMout ) the actual lava flow. Therefore, the higher the FI CUMin and the lower the FI CUMout values are, the more accurate the simulation is (Favalli et al., 2009;Mossoux et al., 2016).

Morphometric modelling
Morphometric parameters and their derivatives of the volcanic cone and lava flow are shown in Tables 1 and 2, respectively. The accuracy of erupted volumes obtained from GIS-based methods is compared with those from classical formulae-derived methods.

Cone
The volcanic cone emerged over a gently N-dipping terrain with a mean slope of 13 ± 6°. The cone displays a horseshoe-shaped crater opening to the NE. Crater rim morphology is influenced by the original relief and wind predominance and shows an asymmetric profile defined by the more significant development of its NW flank, where it reaches a maximum height of 500 m. The eccentricity of the cone (0.7) is smaller than the crater (0.9); therefore, the cone shows a more circular shape in plain view. The mean cone slope is 28 ± 9°, and its internal structure presents a sequence of coarse-grained beds, ranging from welded scoria and volcanic bombs to lapilli. Considering the morphology of this cone and its pyroclastic materials, it can be categorised as type 2 according to Martin and Németh (2006) classification. Therefore, the cone was generated by strombolian eruptive mechanisms (see Figure 3a).
Cone volumes (V DRE ) derived from morphometric modelling and traditional geometric formulae procedures are 563,096 m 3 and 4,851,199 m 3 , respectively, being the latter ∼8.6 times higher. The volume surplus given by the classical formulae-derived procedures is attributed to the poor intrinsic accuracy of this method, which must not be used in volcanic modelling since it could infer volcanic hazard overestimations (Rodriguez-Gonzalez et al., 2010).

Lava flow
The lava flow emission occurs at the crater's cone base (see Figure 3a), from where it is channelled along 2000m following the V-shaped ravine named 'Barranco de Las Martas' (Figure 3b) before reaching a palaeocliff border. Once it is reached, it forms cascades that end up coalescing at the cliff base ( Figure 3c). Then, lava spreads over a narrow abrasion platform and reaches the sea modifying the coastline (Figure  3d-e). Close to the emission point, morphologies of 'a'ā-type lava with levees are present, whereas pahoehoe arch-forming advancement structures are observed through the small lava platform. Lava flow length is 2268 m, including the submarine section, and its width varies from 17 m at the cliff base to 108 m at the end of the Barranco de Las Martas minor incision zone, with a mean value of 54 ± 23 m. Its thickness varies significantly, from 2 m at the cliff border to 21 m at the Barranco de Las Martas major incision zone, where erratic blocks occur, with a  mean value of 5 ± 4 m. Fluvial erosion only affects the N flank since the lava channel itself and the slopes through which it flows facilitate water transport. However, coastal erosion significantly impacts the lava advancement front, causing it to recede several metres, leaving isolated rocks as evidence. Lava flow volumes (V DRE ) obtained from morphometric modelling and traditional geometric formulae procedures are 316,922 m 3 and 459,270 m 3 , respectively, the latter ∼1.4 times higher. Overestimation is also related to the poor intrinsic accuracy of formulae-derived calculations. Thus, GIS-based methods provide more accurate volume estimations either for the volcanic cone or the lava flow.

Lava flow simulations
Calculated and constant standard parameters used as input data for lava flow simulations are presented in Table 3. Q-LavHA does not consider lava flow-seawater interaction (which slows down the lava flow velocity) and the bathymetric information of the DEM, a handicap for correctly simulating lava flows that extend beyond the coastline. However, assuming that the flooding risk of the Montaña de Aguarijo lava flow offshore is null, the DEM used for the simulations in this work is cut and only considers the subaerial path of the flow.

Probabilistic maximum length constraint (L max )
The probabilistic constraint L max depends exclusively on the pre-eruption topography, the lava flow's maximum length, and its minimum and average thicknesses (Table 3). The simulation derived from this approach (see Main Map) is similar to the actual lava flow and shows that the major inundation probability is focused within the central area of the channel. The main zones where the morphology of the simulated lava flow overestimates the real one are also those with a lower inundation probability (see Main Map): (1) near the emission centre, where the simulation runs through a narrow channel semi-parallel to the main flow, (2) at the simulation midpoint, where the NE and SW areas adjacent to the actual lava flow are inundated, and (3) close to the paleocliff termination zone, where the simulation floods the western sectors of the actual lava flow. Contrarily, the underestimated areas are scarce, reduced and mainly located ∼500 m far from the emission centre (SW sector adjacent to the main channel) and in the lava flow termination (E sector adjacent to the main channel). According to the calculated FIs, the overlap (TP), overestimation (FP) and underestimation (FN) percentages of the simulated lava flow are 43.6%, 48.5% and 7.9%, respectively. The composite score (CS) reveals a 26.5% match between the simulated and real lava flow, while the cumulative flood probabilities (CUM) for the pixels located inside (FI CUMin ) and outside (FI CUMout ) the actual flow are 94.1% and 5.9%, respectively. These percentages reveal that the pixels with the highest flood probability match with the actual lava flow and thus that the simulation is correctly adjusted to reality since the input parameters are close to the actual eruption features (Favalli et al., 2009;Mossoux et al., 2016).

Deterministic FLOWGO constraint (L flowgo )
The deterministic model L flowgo considers a set of physical and rheological parameters that determine the lava flow's speed, trajectory, and evolution, in addition to the morphometric parameters already considered in the L max approach. The simulation obtained by this method (see Main Map) shows significant similarities with L max , being the highest flood probability areas located in the central part of the channel flow.
The FIs show a degree of overlap (TP) between the simulated and real lava flows of 42.2%, with 53.7% of pixels overestimated (FP) and 4.1% underestimated (FN). According to the CS, the coincidence between the simulated and real lava flow is 26.5%, with cumulative flood probabilities of 93.6% for the overlapping pixels (FI CUMin ) and 6.4% for the non-overlapping pixels (FI CUMout ). Thus, the L flowgo model, which considers a greater number of parameters than L max , also offers accurate results of the lava flow simulation.
The physical and rheological parameters considered in L flowgo are close to the actual characteristics of the Montaña de Aguarijo lava flow. Therefore, these values could serve as a precedent for the simulation of other lava flows in El Hierro or other intraplate volcanic islands with similar eruptive dynamics since they can be considered the parameters that will control the development of future effusive eruptions.

Conclusions
The compilation of maps of Montaña de Aguarijo's eruption (El Hierro, Canary Islands) illustrates the accuracy of the proposed methodology. Field mapping and generation of 2D and 3D reconstructions allow calculating Montaña de Aguarijo morphometric parameters used as input data to simulate the lava flow inundation probability. The comparison between the resulting probabilistic and deterministic maps shows a high level of correspondence, demonstrating the methodology's usefulness in assessing and forecasting the areas potentially affected by future eruptions on the island.

Software
Georeferencing and digitisation of the eruption were performed initially using GeoPads equipped with the Digital Field Mapping System FieldMove to trace over and vectorise the contact lines. Data processing and development necessary to elaborate the map were performed using TNTGis 2019 and QGIS 3.10 and 3.16 with the plugin Q-LavHA. The 3D images on the map were produced using the TNTGis 2019. The GeoMap application was used for bathymetric maps and Microsoft Excel 2016 to calculate morphometric and simulation lava flow parameters. The topographic profiles were produced using the profile tool QGIS 3.16. The final design was carried out using Adobe Illustrator 2020.

Data availability statement
The authors confirm that the data supporting the findings of this study are available within the article. Acknowledgements Financial support was provided by Project LAJIAL (ref.

Geolocation information
PGC2018-101027-B-I00, MCIU/AEI/FEDER, EU). This study was carried out in the framework of the Research Consolidated Groups GEOVOL (Canary Islands Government, ULPGC) and GEOPAM (Generalitat de Catalunya, 2017 SGR 1494). We thank the collaboration of S. Díaz Rodríguez with the UAV images. We are grateful to the editor M. J. Smith for the effective handling of the manuscript and Makram Murad-al-shaikh, J. Galve and G. Groppelli for providing insightful comments on a preliminary version of this work.