Emplacement conditions of the 1256 AD Al-Madinah lava flow field in Harrat Rahat, Kingdom of Saudi Arabia — Insights from surface morphology and lava flow simulations

Journal of Volcanology and Geothermal Research 309 (2016) 14–30

Contents lists available at ScienceDirect

Journal of Volcanology and Geothermal Research journal homepage: www.elsevier.com/locate/jvolgeores

Emplacement conditions of the 1256 AD Al-Madinah lava flow field in Harrat Rahat, Kingdom of Saudi Arabia — Insights from surface morphology and lava flow simulations Gábor Kereszturi a,⁎, Károly Németh a, Mohammed R Moufti b, Annalisa Cappello c, Hugo Murcia d,e, Gaetana Ganci c, Ciro Del Negro c, Jonathan Procter a, Hani Mahmoud Ali Zahran f a

Volcanic Risk Solutions, Institute of Agriculture and Environment, Massey University, Private Bag 11 222, Palmerston North, New Zealand Geological Hazard Research Unit, King Abdulaziz University, Jeddah, Saudi Arabia c Istituto Nazionale di Geofisica e Vulcanologia, Sezione di Catania, Piazza Roma 2, 95125 Catania, Italy d School of Environment, The University of Auckland, Private Bag 92019, Auckland Mail Centre 1142 Auckland, New Zealand e Departamento de Ciencias Geológicas — Instituto de Investigaciones en Estratigrafía, Universidad de Caldas, Manizales, Colombia f National Center for Earthquakes and Volcanoes, Saudi Geological Survey, Saudi Arabia b

a r t i c l e

i n f o

Article history: Received 1 July 2015 Accepted 2 November 2015 Available online 10 November 2015 Keywords: Scoria cone Basalt Lava flow Lava channel Effusive curve Volume

a b s t r a c t Lava flow hazard modelling requires detailed geological mapping, and a good understanding of emplacement settings and the processes involved in the formation of lava flows. Harrat Rahat, Kingdom of Saudi Arabia, is a large volcanic field, comprising about 1000 predominantly small-volume volcanoes most of which have emitted lava flows of various lengths. A few eruptions took place in this area during the Holocene, and they were located in the northern extreme of the Harrat Rahat, a close proximity to critical infrastructure and population living in Al-Madinah City. In the present study, we combined field work, high resolution digital topography and morphometric analysis to infer the emplacement history of the last historical event in the region represented by the 1256 AD Al-Madinah lava flow field. These data were also used to simulate 1256 AD-type lava flows in the Harrat Rahat by the MAGFLOW lava flow emplacement model, which is able to relate the flow evolution to eruption conditions. The 1256 AD lava flow field extent was mapped at a scale of 1:1000 from a high resolution (0.5 m) Light Detection And Ranging (LiDAR) Digital Terrain Model (DTM), aerial photos with field support. The bulk volume of the lava flow field was estimated at 0.4 km3, while the source volume represented by seven scoria cone was estimated at 0.023 km3. The lava flow covered an area of 60 km2 and reached a maximum length of 23.4 km. The lava flow field comprises about 20.9% of pāhoehoe, 73.8% of 'a'ā, and 5.3% of late-stage outbreaks. Our field observation, also suggests that the lava flows of the Harrat Rahat region are mainly core-dominated and that they formed large lava flow fields by amalgamation of many single channels. These channels mitigated downslope by topography-lava flow and channel–channel interactions, highlighting this typical process that needs to be considered in the volcanic hazard assessment in the region. A series of numerical lava flow simulations was carried out using a range of water content (0.1–1 wt.%), solidification temperature (800–600 °C) and effusion curves (simple and complex curves). These simulations revealed that the MAGFLOW code is sensitive to the changes of water content of the erupting lava magma, while it is less sensitive to solidification temperature and the changes of the shape of effusion curve. The advance rate of the simulated lava flows changed from 0.01 to 0.22 km/h. Using data and observations from the youngest volcanic event of the Harrat Rahat as input parameters to MAGFLOW code, it is possible to provide quantitative limits on this type of hazard. © 2015 Elsevier B.V. All rights reserved.

1. Introduction A compelling reason for studying the evolution of lava flows is to construct predictive models of their behaviour, motivated by the ⁎ Corresponding author at: New Zealand Centre for Precision Agriculture, Institute of Agriculture and Environment, Massey University, Private Bag 11 222, Palmerston North, New Zealand. Tel.: +64 21 1652515. E-mail address: [email protected] (G. Kereszturi).

http://dx.doi.org/10.1016/j.jvolgeores.2015.11.002 0377-0273/© 2015 Elsevier B.V. All rights reserved.

necessity of assessing the immediate hazards posed to people and property by advancing lava flows. Crucial for hazard assessment are estimation of flow paths, flow advance velocities, and final flow lengths that might be expected in future eruptions (Rowland et al., 2005; Wright et al., 2008; Del Negro et al., 2013; Cappello et al., 2015a; Harris and Rowland, 2015). These parameters are easily determined for observed eruptions; however, they are more difficult to infer from solidified flows that represent most of the eruptive history of volcanoes worldwide. For this reason, new approaches have been developed for

G. Kereszturi et al. / Journal of Volcanology and Geothermal Research 309 (2016) 14–30

inferring emplacement conditions from older, solidified lava flows, focusing particularly on flow geomorphology, morphometric measurements of flow units and numerical simulations of lava flows paths (Rowland and Walker, 1987; Rossi, 1997; Soule et al., 2004; Tarquini et al., 2012; Kereszturi et al., 2014; Murcia et al., 2014; Cappello et al., 2015b; Dietterich et al., 2015). The dynamics and emplacement of lava flows are controlled by different parameters, such as viscosity (determined by temperature, chemical composition, gas content, and crystallinity), effusion rate at the vent(s), eruption duration, and the ground topography over which the lava flows (Rowland and Walker, 1990; Griffiths and Fink, 1993; Pinkerton and Wilson, 1994; Cashman et al., 1999; Kerr et al., 2006; Harris and Rowland, 2009). The interplay between such parameters imprints signatures on the lava flows that may vary on multiple spatial scales, and are defined by existent lava surface morphologies, structures and textures. Lava structures include lava rises, lava tubes, tumuli, breakout flows, among many others (Rowland and Walker, 1987; Rossi, 1997; Soule et al., 2004; Guilbaud et al., 2005; Murcia et al., 2014), occurring at a kilometre to metre scale. Lava textures refer to smooth, spiny, blistered, clinker, agglutinated surfaces as well as vesicle distribution (Rossi, 1997; Guilbaud et al., 2005; Murcia et al., 2014), occurring at a much finer scale of metres to millimetres. As a whole, these lava structure and texture characterize the surface morphology that is related to the classical definition of a lava flow after solidification, such as pāhoehoe, 'a'ā and blocky (Macdonald, 1953; Swanson, 1973; Rowland and Walker, 1990; Rossi, 1997; Self et al., 1998). The distribution of flow surface morphologies varies in both space and time, is

15

directly linked to changing flow-emplacement conditions, and results from changes of many internal and external influences down-flow. Despite this, however, quantifiable surface characteristics, defining lava types, are limited by the difficulties inherent in collecting accurate field data on flows with rough topography and large spatial extents. The increasing availability of Light Detection And Ranging (LiDAR) systems allows high-resolution (b 1 m horizontal) topographic data to be obtained, thereby offering opportunities to better understand geomorphic processes from topographic signatures and their hazards (Favalli et al., 2010; Deardorff and Cashman, 2012; Kereszturi et al., 2012b; Cashman et al., 2013). Consequently, the combination of highresolution topography and morphometric analysis can hold important clues to emplacement conditions of already solidified lava flows, when direct measurements of active flows are not possible (e.g. Griffiths, 2000; Soule et al., 2004; Cashman et al., 2013). Significant aspects of lava flow emplacement can also be explored using numerical models that simulate flowing lava. Great advances have been made in understanding the physical processes that control lava flow dynamics and emplacement, as well as the complex feedbacks between cooling, crystallization and rheology that govern those dynamics (e.g. Griffiths, 2000). As a result, over the past years various numerical codes have been developed to predict lava flows footprint and emplacement dynamics (Felpeto et al., 2001; Harris and Rowland, 2001; Crisci et al., 2004; Favalli et al., 2005; Hidaka et al., 2005; Vicari et al., 2007; Connor et al., 2012). However, existing physics-based models cannot yet consider the entire complexity of lava properties, since they have all necessarily involved simplifications of the thermo-

Fig. 1. Overview satellite image (LANDSAT ETM+ 7; R = band 3, G = band 2, B = band 1) of the northern Harrat Rahat, showing the outline of the 1256 AD lava flow field. The red arrows show the lava flows of the Five-fingers eruption. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

16

G. Kereszturi et al. / Journal of Volcanology and Geothermal Research 309 (2016) 14–30

Fig. 2. Effusion sequence of lava flow field development (modified from Camp et al., 1987).

rheological lava properties and relations. This raises the need to evaluate the accuracy and robustness of model results. Also important for hazard assessment is the recognition that, for a lava of given composition, the effusion rate exerts a strong control on the maximum length a volume of lava can attain before it solidifies (e.g. Walker, 1973; Calvari and Pinkerton, 1998). As such, simulations that include the way in which effusion rate changes during an eruption and how this affects lava spreading, are of great interest, as effusion rates can be highly variable over time (Harris et al., 2000; Herault et al., 2009). Considering the abovementioned, MAGFLOW is a useful code that was developed starting from a physical model that includes the rheology of lava and its thermal evolution, including the time-varying effusion

rate as input parameter for lava flow simulations (Del Negro et al., 2008). For example in the MAGFLOW code, the formulations governing the thermal and rheological evolution of a down-flowing lava have a sound physical basis. It was also designed to be adaptive, so that thermal and rheological relationships can be adjusted on the basis of the compositional, rheological or thermal conditions appropriate to the lava flow investigated. By changing input parameters within known and reasonable boundaries and to fit lava flow simulation outputs to field data, it is possible to achieve accurate MAGFLOW simulations of already solidified lava flows. Moreover, numerical simulations driven by timevarying effusion rates not only can provide important insights into lava flow emplacement processes, but they can also be used to estimate

G. Kereszturi et al. / Journal of Volcanology and Geothermal Research 309 (2016) 14–30

17

Fig. 3. High resolution detailed view of the hillshade image derived from the LiDAR DTM with a field photograph of the lava types, such as pāhoehoe (A–B) and 'a'ā (C–D), as well as lava flow structures, such as late stage outbreaks (E–F). The field photo shown in (F) is taken elsewhere than the hillshade image in (E).

the area of the flow field, the lava thickness, and the time required for a flow to reach a particular point and to model topographical changes due to effusive volcanic activity. In this investigation, the emplacement conditions of basaltic lava flows during the 1256 AD eruption at the Harrat Rahat (or Rahat volcanic field), Kingdom of Saudi Arabia, are parameterized through analysing and interpreting high-resolution topographic data. This

information is consequently used as input data to the MAGFLOW code (Del Negro et al., 2008) in order to establish a realistic emplacement conditions that can be used to simulate similar eruption scenarios for lava flow hazard mapping. The 1256 AD event is the youngest manifestation of volcanism at Harrat Rahat. It formed an extensive lava flow field sourced from a fissure in the northern part of the volcanic field in the proximity of Al-Madinah City (Camp et al., 1987; Moufti et al.,

18

G. Kereszturi et al. / Journal of Volcanology and Geothermal Research 309 (2016) 14–30

Fig. 4. (A) Overview hillshade map based on the LiDAR DTM of the source scoria cones of the 1256 AD lava flow field. (B) Slope angle map of the scoria cones. The triangles show the locations used to create the first-order trend surface to model the pre-eruptive surface underneath the scoria cones. (C) Topographic profiles of scoria cones along lines shown in part A.

2012, 2013; Moufti and Németh, 2013; Murcia et al., 2014). Al-Madinah City has now expanded onto these 1256 AD lavas, increasing exposure to any potential renewed volcanism. Most of the Holocene eruptions in the Harrat Rahat produced very similar eruption styles (e.g. lavafountaining to Strombolian type eruptions), dominated by large basaltic lava flows (Moufti and Németh, 2013; Murcia et al., 2014). Therefore, understanding the emplacement conditions of the 1256 AD lava flows and their numerical simulations over space and time may aid assessment of lava flow hazards in this volcanic field. 2. Geological setting and methods 2.1. Description of volcanic activity Harrat Rahat is an extensive intra-plate volcanic field with a dimension of 75 × 310 km bracketed between 24°31′8″N, 39°42′15″E and 21°51′39″N, 40°22′49″E. It formed as part of the Cenozoic intra-plate volcanism in the Arabian Peninsula that was dominated by alkali basalt and hawaiites (Camp and Roobol, 1989; Moufti et al., 2012, 2013). This volcanic field is situated on the Precambrian basement rock, which

comprises plutonic, volcanic and sedimentary rocks (Pellaton, 1981; Stoeser and Camp, 1985). The earliest eruptive product of Harrat Rahat dates back to 10 Ma ago, while the latest eruption is known to have occurred during the 13th century (Camp and Roobol, 1989; Ambraseys et al., 1994). Structurally the region is dominated by N–S and NW–SE normal fault systems related to a failed rift system along the Red Sea (Camp and Roobol, 1989). The Harrat Rahat region erupted in three distinct phases of basaltic units based on K–Ar, Ar–Ar and stratigraphic constraints (Camp and Roobol, 1989; Moufti et al., 2012, 2013; Murcia et al., under review): Shawahit (~ 10–2.5 Ma), Hammah (~ 2.5–1.7 Ma), and Madinah (~ 1.7 Ma—Recent). The latter can be subdivided into a further age-based grouping using the erosional state of the lava flows and volcanic edifices, radiometric and 3He surface exposure ages, as well as stratigraphic and archaeological relationships (Camp and Roobol, 1991; Moufti et al., 2013; Murcia et al., 2015). The 1256 AD event belongs to the youngest age group, along with the neighbouring Five-fingers lava flow field (Fig. 1). The 1256 AD eruption site is located about 20 km south of Al-Madinah City (Fig. 1). The eruption produced a series of scoria cones with an extended (≤10 km from source) tephra blanket (e.g. Kawabata et al.,

G. Kereszturi et al. / Journal of Volcanology and Geothermal Research 309 (2016) 14–30

19

Fig. 5. Closer overview map of the 1256 AD lava flow field, with the location of the two topographic profiles. (B) Topographic profiles showing multiple lava channels and drain-out features.

2015), as well as an extensive lava flow field with a maximum length of 23 km (Figs. 1 and 2) (Camp et al., 1987), which were both deposited on top of previously emplaced lava flows and older tephras. The eruption occurred along a ca. 2-km-long fissure that was initiated by strong earthquakes with a magnitude of up to 4.4 in June 1256 AD (El-Masry et al., 2013), and later lava-fountaining and Strombolian activity formed a chain of scoria cones along the fissure. The chronology of 1256 AD seismic and volcanic events has been reconstructed using historical accounts, based on eyewitnesses and descriptions of local historians (for a review see El-Masry et al., 2013). The extrusion of lava during the 1256 AD eruption occurred in six pulses (pulse 1 = earliest, and pulse 6 = latest events), defined by stratigraphic position of lava flows in Camp et al. (1987). Each pulse resulted in a complex lava flow field and with complex surface morphologies (Camp et al., 1987; Murcia et al., 2014). According to Camp et al. (1987), these pulses of lava extrusion events were fed from three geochemically distinct magmas that were defined as high-K2O, low-K2O end-members, and a hybrid type that is an intermediate type between the two endmembers. It is noteworthy that Murcia et al. (under review) proposes a different interpretation where the magmas are not distinct and rather represent the evolution of a single one. The pulse 1 fed by both highK2O, low-K2O end-members magmas formed two channels dominated by pāhoehoe surfaces that travelled a maximum distance of 13 km northwards of the source vent (Fig. 2). This pulse formed two lava branches (Fig. 2). The pulse 2, which is about 70% of the total emitted volumes, formed by the three high-K2O, low-K2O end-members, and hybrid type magmas (Camp et al., 1987). The pulse products are superimposed on top of the deposits of the previous flow. This lava flow reached an area

that is about 12 km northward from the source fissure, and then it split to two branches (Fig. 2). One branch reached a maximum distance of 23 km north of the source vents. The pulse 3 formed two lava channels in which one was directed to N (length = 9.5 km), following the lava flows of pulse 1 and 2, while one was directed to NW (length = ca. 5 km), forming the narrow branch 1 by juxtaposition of the deposits from pulse 1 (Fig. 2). The latest lava flow from pulse 4 to 6 emplaced about 5% of the total emitted volume. These late-stage flows contributed little to the overall length and geometry of the lava flow field (Camp et al., 1987). Pulse 4 reached a length of 4.5 km, while pulse 5 and 6 formed maximum length of ca. 3.5 km and 2 km, respectively (Fig. 2). The morphotypes of these lavas are pāhoehoe dominated (Camp et al., 1987). In such proximal areas, agglutinated rafts of scoria, originated from the source scoria cones, are common on lava flow (Camp et al., 1987).

2.2. Geological mapping The outline of the scoria cones and lava flow from the 1256 AD eruption was mapped out at a scale of 1:1000 in a Geographical Information System (GIS), using field notes and field data, high resolution orthophotos (0.6 m), Landsat 7 ETM+ satellite imagery (15–30 m resolution), a high resolution (0.5 m) Light Detection and Ranging (LiDAR) Digital Terrain Model (DTM) and its derivatives, such as hillshade and slope angle maps, and previous work by Camp et al. (1987). This mapping was also supported by field-based verification of ambiguous locations during multiple field campaigns. Individual lava

20

G. Kereszturi et al. / Journal of Volcanology and Geothermal Research 309 (2016) 14–30

Fig. 6. Branches (A) and flow perpendicular transects of the 1256 AD lava flow (B).

flow units and lava flow structures, such as channels and lava lobe, were digitized manually. In this study, a new mapping technique was used to establish the proportion and spatial locations of dominant lava types as a function of distances from their sources. Using the digitized polygon of the lava flow field, 1500 randomly distributed points with a minimum allowed spacing of 75 m was used. At these point locations, the lava types, such as pāhoehoe (Fig. 3A–B) and 'a'ā (Fig. 3C–D), as well as lava flow structures, such as late stage outbreaks (Fig. 3E–F), were assessed manually using a high-resolution hillshade image derived from the LiDAR DTM.

(2012a). This approach was preferred due to the potential burial effects of lava flow activity (e.g. Favalli et al., 2009) during and after scoria cone emplacement. Modelling the pre-eruptive surface underneath the scoria cone is reasonable due to the overall flat topography (slope angle is ≤1°). The pre-eruptive surface was defined as a first-order trend surface created from digitized spot heights (n = 300) outside of the lava flow units (Fig. 4). The volume is defined as a 3D space enclosed by this modelled pre-eruptive surface and the LiDAR DTM within the area of mapped scoria cones (Fig. 4).

2.3. Modelling subsurface geometry and eruptive volumes The eruptive volume calculations for the scoria cones were carried out, using the inclined plane approach, similar to Kereszturi et al. Table 1 MAGFLOW input parameters for the 1256 AD lava flow simulations. Input parameters

Symbol

Unit

Value used

Density Initial temperature Solidification temperature Emissivity Specific heat H2O content Bulk volume Effusion curve Duration Source

ρ Ti Ts ε cv – v – t –

kg/m3 C C – J/kg/K wt.% km3 – days –

2700 1150 800–600 0.9 1150 0.1–1 0.4 Peaked/constant 52 Single/multiple

Fig. 7. Graph showing the two effusion curves used in the MAGFLOW simulations.

G. Kereszturi et al. / Journal of Volcanology and Geothermal Research 309 (2016) 14–30 Table 2 Volume estimates for the 1256 AD scoria cones. Parameters

Unit

Cones 6–7

Cone 5

Cone 4

Cone 3

Cone 2

Cone 1

Surface area Bulk volume

106 m2 106 m3

0.41 14.80

0.12 3.69

0.17 4.70

0.02 0.18

0.01 0.11

0.02 0.16

Table 3 Geometry of the 1256 AD lava flow field. flow segment

Maximum length

Branch 1 Branch 2 Branch 3 Total

Area 2

m

km

23.4 5.8 5.4 23.4

ca. 48.7 ca. 2.9 ca. 8.2 59.8

Bulk volume 106 m3 329.8 9.8 51.9 392

The lava flow field emplaced on an irregular topography, dominated by pre-1256 AD lava flows and older volcanic edifices (e.g. Fig. 5A and B). In order to rebuild this complex pre-eruptive morphology, a Triangulated Irregular Network (TIN) model was used, created from spot heights spaced 1 m apart and located 2 m from the present flow margin. Elevations at these spot heights were extracted from the LiDAR DTM and a TIN model was constructed using Delaunay criteria. The general pre-eruptive terrain slope was also established using a first-order trend surface fitted to 1000 randomly generated points located along the digitized lava flow margin. 2.4. Lava flow geomorphology and braiding index The 1256 AD lava flow field can be subdivided into three branches (Fig. 6A–B) for the ease of computation and data-handling. The lava flow field boundary was used to extract centre line of each branch. The vertices along the polyline were densified with 1 m point-spacing and were then converted to points. Those points along the flow margin were subsequently used to construct Thiessen polygons. These polygons divide the available space between input points, providing the centre line of the combined flow unit (Fig. 6A). The centre line was then smoothed using a Polynomial Approximation with Exponential Kernel algorithm with a kernel size of 1000 m. This ensures that the extracted centre line is free of artificial sharp bends (Fig. 6A). Along this smoothed centre line, perpendicular lines were constructed at 50 m spacing (Fig. 6B). Along these transects, flow geometry (e.g. total flow width) and average and maximum thicknesses were calculated. The thickness values were extracted from the thickness map, produced from the extraction of the pre-eruptive surface from the present LiDAR DTM, using zonal

21

statistics along each perpendicular transect. Slope angles were measured in a 100-m-wide strip outside of the lava flow. Average and median slope angles of the pre-eruptive terrain were calculated for the entire LiDAR covered area, and then later masked to a 100-m-wide strip outside of the lava flow in order to avoid edge effects. The slope angle was calculated with a finite difference algorithm using an unweighted, linear Prewitt filter (e.g. Prewitt, 1970). The braiding index was calculated in a similar way to Egozi and Ashmore (2008) and Dietterich and Cashman (2014). The braiding index shows the number of branches, defined as a central line of a lava flow channel, intersected at the flow transect lines. Thus, the larger braiding index indicates that a larger number of potential contributing channels developed at a given down-flow segment of a lava flow field. The subdivision of lava flow into three branches helped avoiding artificial increase of braiding index due to bending of the lava flow and complex outline of the 1256 AD flow field. 2.5. Numerical modelling with MAGFLOW The MAGFLOW code for lava flow simulations is based on the cellular automata method for 2D computation systems, using a spatial partition into cells whose state evolves according to that of their neighbours and a transition function that describes the exchanges of lava and heat (Del Negro et al., 2008). The MAGFLOW code evaluates the lava exchange between neighbouring cells from the steady state solution of Navier–Stokes equations for the motion of a Bingham fluid on an inclined plane subject to pressure force, in which the conservation of mass is guaranteed both locally and globally (Del Negro et al., 2008; Cappello et al., 2015a). MAGFLOW code requires the following input parameters: the rheology (eruption and solidification temperatures, lava density, emissivity, specific heat, water content), effusive curve (i.e. effusion rate estimates over a given time period), source location, and a high-resolution DTM of the terrain over which the lava flows. A list of input parameters used in the MAGFLOW numerical simulations is tabulated in Table 1. The rheological components were constrained using geochemical properties of the 1256 AD lava flow field after Moufti et al. (2012) and Murcia et al. (under review). This lava flow has a tight range in whole-rock SiO2 content (45.60–46.84 wt.%) with a wider range in MgO content (9.29 and 5.50 wt.%; Mg# 63–46), which reflects fractionation of clinopyroxene and olivine en-route to the surface. The lava flow was fed by alkaline basaltic magma. Field constraints suggest that the eruption progressed from less to more evolved products, with the last part of the eruption being the most explosive (Camp et al., 1987; Murcia et al., under review). The major mineral phases are olivine (Ol) and plagioclase (Pl) in the lava rocks and scoriaceous deposits, while minor mineral phases are clinopyroxene and titanomagnetite (Camp

Table 4 Comparison of 1256 AD lava flows to other well-known modelled and observed lava flow geometries and physical characteristics. Parameters

Parameters Age

Duration

Maximum length

Area

Bulk volume

Flux Average

Etna Mt. Roskill Jebel at Tair Three Kings Krafla Collier flow Lonquimay Medinah flow Paricutin Carrizozo Laki

–

Days

km

km2

km3

m3/s

2011 30 ka 2007–08 AD 28 ka 1984 AD 1600 yr BP 1988–90 AD 1256 AD 1943–52 AD 2 ka 1783–84 AD

1.5 21 ca. 115 83 14 255–275 330 52 3293 – 245

4.3 2.5 ca. 2 6.7 9 13.6 10.2 23.4 ca. 5 75 35

1.0 2.0 5.9 7.6 12.1 – – 59.8 25 – ca. 600

0.001 0.015 0.022 0.078 0.134 0.170 0.23 0.392 0.7 4.5 14.7

10 8 2 11 51–91 10–50 8 87 ca. 2–3 3100–11000 8500–8700

Reference Maximum

60 40 – 40 162 38 – 395–730 6 – –

Vicari et al. (2011) Kereszturi et al. (2014) Xu and Jónsson (2014) Kereszturi et al. (2014) Rossi (1997) Deardorff and Cashman (2012) Kerr and Lyman (2007) this study Krauskopf (1948), Fries, 1953 Keszthelyi and Pieri (1993) Thordarson and Self (1993)

22

G. Kereszturi et al. / Journal of Volcanology and Geothermal Research 309 (2016) 14–30

Fig. 8. A series of maps showing: (A) lava flow thicknesses modelled using a TIN-based method, (B) point counting results of lava type/structure, (C) lava flow channels and their centre lines and, (D) braiding index.

et al., 1987). The phenocrystals are present up to 8 vol.% (Ol ≈ Pl), but are commonly ≤2 vol.% (Ol N Pl) for the latter eruptive products. The origin of the magma is estimated to be ~ 100 km deep, in the asthenospheric mantle, with temperatures of ~ 1500 °C and pressures of N3 GPa. The eruption temperature is estimated to be either between

Fig. 9. Histogram of thickness values calculated from the LiDAR DTM with the average, median and model thickness values of the 1256 AD lava flow field.

1130 °C and 1230 °C, based on MgO geothermometers (e.g. Helz and Thornber, 1987), or close to, but no lower than, 1100 °C applying plagioclase thermometry. As an input parameter for MAGFLOW simulation, an intermediate value of 1150 °C was used. Based on the sensitivity analyses of MAGFLOW code, Bilotta et al. (2012) concluded that solidification temperature, water content and effusion curves are the most sensitive input parameters of the lava flow simulations. Hence, in the present study a range of input parameters were used. The solidification temperature is modelled to occur at 800 °C, 700 °C and 600 °C (Table 1). As the extrusion temperature was kept constant, the solidification temperatures above correspond to a total cooling of 350 °C, 450 °C and 550 °C, which is in line with observation of surface temperature and flow advance patterns on Kilauea, Hawai'i (e.g. Ball et al., 2008). Direct measurements of water content were carried out on olivine crystals melt inclusions from the lava flow with Fourier Transform Infrared Spectroscopy (FTIR). This returned a water content range from 0.2 to 1 wt.% (Murcia, H. unpublished data). Therefore, an order of magnitude change is implemented into the MAGFLOW lava flow simulations, by inputting water content values of 0.1, 0.5 and 1 wt.% (Table 1). The simulations above carried out from a single point-like source, however, there is evidence that the effusive activity of the 1256 AD Al-Madinah eruption lasted about 52 days (1248 h) based on historical descriptions (Camp et al., 1987; El-Masry et al., 2013). This long-lasting eruption is rather formed from pulse of lava flows than a single lava flow, evidenced by superimposed lava flow channels and stratigraphic position of lava flow units (e.g. Fig. 2). To account for such complex lava flow field development two types of effusion curves were used: simple and complex (Fig. 7). The “simple” effusion curve has a linear increase with a break in slope at 1/8 of the entire duration of 52 days

G. Kereszturi et al. / Journal of Volcanology and Geothermal Research 309 (2016) 14–30

23

Fig. 10. Variability of flow width (black line), underlying slope angle (blue line) and lava flow thickness (red line) parameters in a down slope direction along the main branch 1 (A) and the two peripheral branch 2 and 3 (B and C). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

(Fig. 7), where after the peak magma flux there is an exponential decay until 1 m3/s, similar to Kereszturi et al. (2014). This one was used to conduct a sensitivity analysis on the effect of solidification temperature and water content. The “complex” effusion curve has similar shape to simple one, but it is portioned into six individual flow units each with a total of 8 days length (Fig. 7), consistent with field-based evidence (Camp et al., 1987). The estimated total bulk volume from the LiDAR DTM was used, and it was proportioned based on Camp et al. (1987). These flow units were simulated sequentially, meaning that after each simulation the input DTM of the subsequent simulation was updated with the effusion

product of the previous flow unit, ensuring lava flow channel–channel interactions to occur during emplacement. Lack of independent data (e.g. geophysical profiles and drill cores) regarding the morphology of the pre-eruptive terrain underneath the 1256 AD lava flow field contributes to the uncertainties of producing an accurate simulation that recreates the emplacement of this flow as a whole. However, the information derived from the 1256 AD flow field can easily be analysed to provide inputs into the MAGFLOW code. These inputs can then be applied to simulate a lava flow anywhere in the Harrat Rahat to be able to give an idea of future lava flow hazards

24

to

G. Kereszturi et al. / Journal of Volcanology and Geothermal Research 309 (2016) 14–30

Al-Madinah City. Our approach to date has been to run a simulation of a future eruption scenario that assumes the source or vent is 2 km south of the original eruption site. The new vent is located as close as possible to the 1256 AD vents but is placed in a location that ensures that past lava flows (N1256 AD) are not a barrier to the advancement of the simulated lava flows. The simulation was carried out on the present topography represented by a resampled LiDAR DTM with a horizontal resolution of 10 m. 3. Results and interpretation 3.1. Eruptive volumes and extrusion rates The total bulk volumes of the source scoria cones are 23 × 106 m3, which is made up from seven scoria cones, each with relatively small bulk eruptive volumes (Table 2). This volume is in the same range as other scoria cones from the Harrat Rahat estimated from 30 m resolution ASTER GDEM (cf. Zaidi and Mukhopadhyay, 2015). The combined volume of the lava flow field is at least two orders of magnitude larger than the scoria cone volumes, which are estimated to be 392 × 106 m3 (Table 3). This lava flow volume is in agreement with previous estimates of 290 × 106 m3 from Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) DTM (Murcia et al., 2014), and ca. 500 × 106 m3 from the area multiplied by average thickness method (Camp et al., 1987). This range of bulk volumes for lava flows often occurs in other volcanic fields, such as Jorullo, Mexican Volcanic Belt (Guilbaud et al., 2011), Southeast Crater Flat, Southwestern Nevada volcanic field, U.S.A. (Valentine and Perry, 2006), Rangitoto, Auckland volcanic field, New Zealand (Kereszturi et al., 2013), and Agár-tető, Bakony–Balaton Highland volcanic field, Hungary (Kereszturi and Németh, 2012). The estimated scoria cone and lava flow bulk eruptive volumes, combined with the emplacement duration of 52 days, result in an overall extrusion rate of ca. 5 m3/s for the scoria cone building phase, and an average effusion rate of 87 m3/s for the period of lava flow emplacement (Table 4). This range of average extrusion rates compares well with other well-studied lava flows from flank eruptions and on monogenetic volcanic fields; the 1256 AD lava flow field represents a medium to large volume eruption scenario in such magmatic systems (Table 4). 3.2. Lava flow morphometry A GIS-based morphometric analysis of the entire lava flow field was carried out, determining the combined flow width, thickness and slope angle of the pre-eruptive terrain (Figs. 8–10). The thickness calculated with the TIN-based basement reconstruction yielded a maximum thickness of 30 m (Figs. 8A and 9). These thicknesses are limited to three well-defined geographical areas (Fig. 8A). The average and median thicknesses of the lava flow are, however, 6.6 m and 6.2 m, respectively (Fig. 9). Along the longest branch 1, there is no general distinctive decrease of flow width and maximum thickness (Fig. 10A). However, locally, there are distinctive changes in maximum thickness values (Fig. 10A). This is interpreted to be due to the effect of pre-eruptive topographic obstacles (e.g. pre-existing scoria cones) and depressions (Fig. 9A). These places caused the lava flow front to halt and subsequently divert and split in to narrower lava channels (Kauahikaua et al., 2002). These interactions of lava flow units with the fine topography (e.g. lava flow edges, topographic depressions, or scoria cones) are very prominent features of the Harrat Rahat's lava flows. These lava flows (not just those from the 1256 AD lava flow field) hence commonly emplace “laterally” (e.g. Fig. 5B), whereby one flow unit flows and interacts with topography as described above; the flow unit comes to rest Fig. 11. Histograms showing the distribution of lava types, such as pāhoehoe (A) and 'a'ā (B), as well as lava flow structures, such as late stage lava outbreaks (C) from proximal to distal locations.

G. Kereszturi et al. / Journal of Volcanology and Geothermal Research 309 (2016) 14–30 Table 5 Input bulk volumes used in the MAGFLOW simulations. The volume proportions are after Camp et al. (1987). Pulse

Simulation parameters V

Proportion 6

3

10 m 0.03 0.278 0.05 0.018 0.019 0.004 0.40

Units Pulse 1 Pulse 2 Pulse 3 Pulse 4 Pulse 5 Pulse 6 Combined

% 7.5 69.5 12.6 4.6 4.7 1.1 100

and then the next one usually flows next to this already immobile flow, causing a lateral spread. Branch 1 has a maximum length of 23.4 km, while branch 2 and 3 are the two left-lateral arms with lengths of 5.8 km and 5.4 km, respectively (Table 3). The channel margins and centre lines (Fig. 8C) were digitized directly from the hillshade image, and were used to establish a braiding index for the 1256 AD lava flow field (Fig. 8D). The maximum braiding index was calculated to be 6 (Fig. 8D). Other studies show much larger braiding indices, such as the 1984 AD Mauna Loa lava flow which had a braiding index of up to 12 (Dietterich and Cashman, 2014). This can be explained by changes in lava viscosity down-flow, as the more viscous the lava becomes the smaller the braiding index (e.g. Dietterich and Cashman, 2014). The lava flow width changes locally along the entire flow field. Although, there is no general trend, which is different to other lava flow, such as Collier Cone lava flow in Oregon, U.S.A. (e.g. Deardorff and Cashman, 2012). The lava flow morphometry changes more clearly in a down-flow direction for branches 2 and 3 than for the branch 1, with a prominent decrease in flow width and maximum thickness (Fig. 10B and C). These two branches are made of a single (branch 3) and two narrow lava channels (branch 2). These morphometry changes are interpreted here as evidence that these marginal branches were shorter-lived than the main branch 1; consequently the original channel morphology could be preserved (i.e. no overprinting of processes). In these two marginal branches, the systematic down-flow maximum thickness and width decreases indicate the dependence of them on the available lava supply. Overall, the lack of any obvious correlation of lava flow width with pre-eruptive slope angle (e.g. Fig. 10A) is interpreted to be due to the sequential development of the

25

flow field. This means that the individual channels emplaced at different times during the course of the eruption, forming a complex lava flow field. It is worthwhile noting that the pre-eruptive surface slope does not change dramatically, but it is characterised by a high frequency fluctuation (Fig. 10A–C). This high frequency fluctuation is due to the finescale (metre-scale) topographical features that the LiDAR DTM could resolve. 3.3. Lava types The lava flow types of the 1256 AD flow field change from proximal pāhoehoe to distal 'a'ā (e.g. Fig. 8B). The area proportions of the lava types and lava outbreak structures were established using random point counting on hillshade images as 20.9% of pāhoehoe, 73.8% of 'a'ā, and 5.3% of late-stage outbreaks (Fig. 8B). The pāhoehoe type occurs very close to the source (Fig. 11A), while the 'a'ā type occurs most often in medial to distal parts of the lava flow field (Fig. 11B). The late stage small-volume lava outbreaks occur everywhere along the margins of the flow (Fig. 11C). This spatial pattern of lava morphotypes can be interpreted in two ways. (1) This shows that the surface evolution started with initial pāhoehoe sheet flows, which transformed to 'a'ā-filled channels within pāhoehoe sheets, and then to channelized 'a'ā in distal regions. This is very similar to the evolution of the Five-fingers lava flows in the Harrat Rahat, and Hawaiian 'a'ā flows (Soule et al., 2004). These changes are also consistent with changes in the fluid dynamics from initial Newtonian to Bingham type fluids behaviour (Pinkerton and Wilson, 1994; Castruccio et al., 2013). This transition in lava type (e.g. from pāhoehoe to 'a'ā) occurs in the spatial scale of 5 km down-flow (Fig. 11). This evolution matches perfectly with the transition (ca. 4 km down-flow) that occurred from pāhoehoe to 'a'ā type lava on the December 1974 AD lava flow of Kīlauea, Hawai'i (Soule et al., 2004). (2) This spatial transition of morphotypes can also indicate a temporal change of the eruption due the fact that late-stage lava flow, dominated by pāhoehoe morphotypes superimposed on top of earlier lava flow products. The products of the late-stage pulses (flows 4–6) travelled maximum distance of 5 km from the source area (Fig. 2). Pāhoehoe and 'a'ā lava types can be defined as crust-dominated and core-dominated flow regimes, respectively (e.g. Castruccio et al., 2013). The dominant proportion of the 1256 AD lava flow field was emplaced as an 'a'ā type, highlighting that this lava-flow field formed through a core-dominated regime. This is consistent with the dominance of lava flow channels that developed during

Table 6 Area of inundation and length of simulated lava flow by MAGFLOW. Source type

Single source

multiple Real

Simulation

S1 S2 S3 S4 S5 S6 S7 S8 S9 F1 F2 F3 F4 F5 F6 Merged Combined

Water content

Solidification temperature

Length

Thickness

Area

Simulated

Measured

Maximum

Mean

Simulated

wt.%

°C

km

km

m

m

km2

0.1 0.1 0.1 0.5 0.5 0.5 1 1 1 0.5 0.5 0.5 0.5 0.5 0.5 0.5 ca. 0.5

800 700 600 800 700 600 800 700 600 800 800 800 800 800 800 800 ca. 800

22.37 22.10 21.84 25.63 26.41 25.42 28.66 27.55 30.72 11.43 26.07 12.25 10.94 8.21 4.47 26.07 –

– – – – – – – – – 13 23 10 5 4 2 23 23.4

24.70 31.60 32.87 29.70 31.08 30.15 31.35 31.90 27.10 17.21 30.16 23.42 15.94 14.08 12.42 – 25.17

10.70 13.30 13.25 10.08 10.18 10.42 8.14 8.93 5.32 6.35 9.08 6.80 5.31 4.85 4.40 – 6.2

25.79 27.42 27.06 32.59 33.65 32.06 34.14 34.53 28.73 4.60 27.01 7.26 3.39 3.87 1.02 38.85 59.80

S9 simulation (in bold and italic) reached the edge of the LiDAR DTM coverage, thus these values are considered as minimum.

26

G. Kereszturi et al. / Journal of Volcanology and Geothermal Research 309 (2016) 14–30

emplacement (e.g. Fig. 8C), and with the qualitative model of the spatial distribution of lava types in the 1256 AD flow field (Camp et al., 1987; Murcia et al., 2014). 3.4. Lava flow simulations The MAGFLOW simulations were carried out based on the results on the lava flow morphometry and volumes (Table 3), as well as data on geochemistry, effusion flow units (Fig. 7 and Table 5), and eruption temperature estimates (Murcia et al., 2014, under review). The lava flow simulation from a single source was tested with a range of input values of solidification temperatures and water contents (Table 6). In total nine simulations have been carried out to test the sensitivity of the results to this input parameters. Based on the range of values the best match between simulation and the actual flow extent and flow length is the case of low water content (e.g. S4 in Table 6). If the water content increased from 0.1 wt.% to 1 wt.%, the run-out distance of the lava flow increased significantly (Fig. 12). On the other hand, the changing the solidification temperature from 800 °C to 600 °C had an poor effect on the lava flow length, and a minor effect on the total area of inundation (Fig. 12). All simulated lava flow reach their maximum run-out-distance after ca. 300 h (b12 days) from the initial extrusion, regardless of water content (Fig. 13). This corresponds to a range of lava flow advance velocities between 0.01 and 0.22 km/h. The most prominent feature of the simulated lava flows, regardless of what solidification temperature and water content was used, is the narrow channels in proximal locations (e.g. Fig. 12). This was interpreted to be due to the solidification processes, thus, once the initial lava channel partially solidifies there is not lateral spreading in proximal settings. This results in an underestimation of overall inundation area in proximal locations. The multiple source simulation returned a lava flow outline that is more compatible with the actual 1256 AD lava flow field outlines (Figs. 14 and 15; Table 6) than the single simulations. The total area of inundation is slightly larger than in the single simulations, this is in part due to the fact that the subsequent simulation allowed us to model channel to channel interaction, leading to larger area of inundation in proximal locations (Figs. 14 and 15). The total area of inundation for the simulation with the complex effusion curve is 38.8 km2, while from a single source (with the same solidification and water content parameter) is 32.6 km2 (Table 6). This area of inundation is still about 65 and 55% of the actual extent of the 1256 AD lava flow field, respectively. The difference in the spatial distribution of lava flows is about 16% and shows that in the this special case the effusion curve has only a moderate effect on the resultant simulated lava flow outline and area of inundation in the Harrat Rahat area. 4. Discussion 4.1. Opportunities and challenges in numerical simulations of lava flow MAGFLOW simulations and mapping data on the lava flow's morphology and lava types support that the 1256 AD flow field attained a length that was in accordance with the magmatic controls, such as magma flux and the rheological evolution (e.g. crust growth, cooling, crystal content, ponding of lava). The 1256 AD lava flow field from the Harrat Rahat shows a complex emplacement history with stages of lava ponding and pulsation of magma supply (e.g. increased thickness area in Fig. 8A). This led to complex emplacement processes influenced by the balance between preeruptive topography (e.g. older lava flow and scoria cones) and effusion rate at the lava flow source. This complex behaviour can only partially Fig. 12. Result of MAGFLOW simulation using a single effusion curve (Fig. 7) with changing solidification temperatures (Tsol) and water contents (W).

G. Kereszturi et al. / Journal of Volcanology and Geothermal Research 309 (2016) 14–30

Fig. 13. Change in maximum length during the first 500 h of simulation for S1 (solid blue curve), S4 (dashed red curve) and S7 simulation settings (dash-dot black curve). Note that they reach the maximum run-out-distance at about 300 h after extrusion. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

be modelled using MAGFLOW, which is able to take into account the way in which the effusion rate changes during eruption (e.g. complex effusion curves; Fig. 7), including the discrete number of pulses each with a relatively small eruptive volume. This approach of simulating lava flows is consistent with the petrological and geochemical model of the 1256 AD eruption which involves multiple batches of ascending magma and a complex temporal evolution of the lava flow field

27

(Camp et al., 1987; Murcia et al., under review). Simulating such complex emplacement patterns (e.g. Fig. 14), however, returned a slightly larger area of inundation (Table 6), highlighting that even simulations using a simple, more generalized, effusion curve could be used for simulating lava flows. This is consistent with the sensitivity analyses of the MAGFLOW code, which showed that simulations based on continuous estimates of effusion rate, even with a moderate error range, yield a better result than simulations based on a fewer, but accurate discharge rate data (Bilotta et al., 2011; Ganci et al., 2012). The initial stage of lava emission in the Harrat Rahat region usually takes place as a pāhoehoe type. Moreover, the lava flows often show emplacement structures that are generated by a pulsatory magma supply and lava pressure, leading to formation of collapse pits, pressure ridges and tumuli fields (c.f. Murcia et al., 2014). Notwithstanding the complexity of these conditions, the simulated lava flow with similar parameters to the historical 1256 AD flows showed that the resultant geometry, morphology, length and emplacement processes (e.g. bifurcation, interaction with the micro-topography) are very similar to the actual 1256 AD flows emplacement processes (Fig. 15). On the other hand, proximal features of the simulated lava flows are poorly represented, resulting in a potential underestimation of volcanic hazard in such settings. Based on the simulation carried out in the present study, medial and distal position of the simulated lava flow, and the maximum run-out distance of flows can be used to constrain lava flow hazard maps in the region. All existing flow emplacement models adopt several simplifications to reproduce lava flow dynamics and cannot yet completely simulate lava flow behaviour, such as lava tube and ephemeral vent development and transitions from pāhoehoe to 'a'ā modes, and they must be utilized with precaution. The effective measure of how well lava flows are understood will govern our ability to simulate their advance by numerical modelling. Each model works well for a specific set of constraints, but there is still a need to develop a more robust, generalized approach that fully links the thermal and rheological evolution of the lava to the flow dynamics. In addition, it is worth highlighting that numerical simulations are critically sensitive to uncertainties in input parameters, such as effusion rate, the changing rheology of molten lava as it cools

Fig. 14. Results of MAGFLOW simulations using a complex, pulse-based effusion curve, shown in Fig. 7.

28

G. Kereszturi et al. / Journal of Volcanology and Geothermal Research 309 (2016) 14–30

Fig. 15. Comparison of simulated lava flow field outline with the 1256 AD lava flow field as a reference (black polygon).

down-flow, and topographic attributes of the terrain (e.g. steepness and morphology); propagation of error through the resulting simulated flow may have a non-negligible effect on the accuracy and reliability of any prediction based on the model (Bilotta et al., 2012). This study showed that a change in solidification temperature results in a minor fluctuation of the maximum lava flow length and the area of inundation, while changing the water content could change the resultant lava flow's area by 34%, and maximum length by 28% (e.g. Table 6). This can be mitigated by direct water content estimates of the eruption products, and/ or simulating reasonable range of input parameters, as it has been carried out in the present study. 4.2. Volcanic hazard implications Assessing the hazard of inundation from lava flows in monogenetic volcanic fields or on composite volcanoes is getting a new focus in volcanology due to the increasing availability of high-resolution DTMs, such as LiDAR, and numerical codes for lava flow simulations (Hidaka et al., 2005; Vicari et al., 2007; Favalli et al., 2010; Kereszturi et al., 2012b; Harris and Rowland, 2015). However, the selection of an appropriate physics-based model of lava flow evolution and emplacement is fundamental for accurate prediction of likely paths of lava flows and their areal extents. Previously, lava flow hazard in the Harrat Rahat was mapped out by the topography-based code of Connor et al. (2012), using user-defined thickness values and volumes (El Difrawy et al., 2013). This approach has some limitations for predicting length of lava flows. In these topography-based simulations, the lava flow

length was attained using randomized lava thickness values between 0.05 and 0.3 m. These low values potentially predict longer lava flows with lengths ≥20 km. Using the lava flow simulations, the time-averaged advance rate of lava flow can be calculated for every 10 days. This shows that a lava flow like in the 1256 AD eruption, advanced with an average speed of ca. 0.1 km/h in the first 10 days. The maximum speed reached ca. 0.2 km/h in the lava flow simulations (e.g. S7–9). Realistically, the initial few days the lava flow with a finite magma supply will reach its maximum length extent, as documented in examples from Hawai'i (e.g. Riker et al., 2009). As a result, the MAGFLOW code can be used as a tool to estimate lead-times before lava flow invasion occurs at a critical infrastructure within the Harrat Rahat, increasing pre-eruption preparedness. 5. Conclusion This study focused on the eruption mechanisms and distribution of the 1256 AD lava flows in order to provide quantitative limits on this hazard for the Harrat Rahat. Here, field work combined with highresolution digital topography, morphometric analysis and numerical simulations were used to infer the emplacement history of the 1256 AD Al-Madinah lava flow field. Our main conclusions are summarized in the following points. • Based on the numerical simulation results using the 1256 AD flows' parameters, we can conclude that the MAGFLOW code is able to model rheological changes due to cooling during lava flow

G. Kereszturi et al. / Journal of Volcanology and Geothermal Research 309 (2016) 14–30

•

•

•

•

emplacement. This can successfully predict the run-out distance of lava flows; however, for this code to be efficient requires many input parameters, such as solidification temperature, density, emissivity, specific heat, water content, total volume, effusive curve, source location, and pre-eruptive topography. Cooling and topographic influences on lava flow emplacement processes are very important in the Harrat Rahat region and need to be incorporated appropriately before any volcanic hazard assessment. Pre-eruption evaluation of lava flow hazard in this region is hampered by ponding and the subsequent extrusion cycle of lava flows. This might be mitigated by syn-eruptive topographic surveys, by aerial photogrammetry and LiDAR surveys for a near real-time surveillance of lava flow emplacement. The setting of the 1256 AD lava flow field can be used to simulate lava flows of this volumetric size. This volumetric size might be an extreme in the Harrat Rahat region, thus the systematic parameterisation of lava flows of smaller volumetric sizes, such as the Five-fingers lava flow fields (Fig. 1), must be included in lava flow hazard assessment of the Harrat Rahat. This can expand the variety of eruption scenarios that can serve as a basis of a robust hazard estimate of the Harrat Rahat region. LiDAR-generated high-resolution DTMs provide a good means of analysing details of lava flow surfaces at a km to m scales. Due to the available details, flow areas and thicknesses can be delimited from the DTMs, improving estimates of flow volumes of solidified lava flows. Airborne LiDAR can also be used as a rapid mapping tool, providing support data for lava flow simulations, such as accurate eruptive volume and magma supply estimates using multiple-temporal airborne surveys during a volcanic crisis. This option and its implementation to numerical simulation in the Harrat Rahat region require further investigations.

Acknowledgement GK is thankful for the logistic support from the Geohazard Research Center of King Abdulaziz University, Jeddah, Kingdom of Saudi Arabia. The MAGFLOW model was developed within the framework of TechnoLab, the Laboratory for Technological Advance in Volcano Geophysics, at INGV in Catania, Italy. The authors are grateful to MarieNoelle Guilbaud and an anonymous reviewer for the valuable comments and suggestions. References Ambraseys, N.N., Melville, C.P., Adams, R.D., 1994. The Seismicity of Egypt, Arabia and the Red Sea. Cambridge University Press. Ball, M., Pinkerton, H., Harris, A.J.L., 2008. Surface cooling, advection and the development of different surface textures on active lavas on Kilauea, Hawai'i. J. Volcanol. Geotherm. Res. 173 (1–2), 148–156. Bilotta, G., Rustico, E., Hérault, A., Vicari, A., Russo, G., Del Negro, C., Gallo, G., 2011. Porting and optimizing MAGFLOW on CUDA. Ann. Geophys. 54 (5), 580–591. Bilotta, G., Cappello, A., Hérault, A., Vicari, A., Russo, G., Del Negro, C., 2012. Sensitivity analysis of the MAGFLOW cellular automaton model for lava flow simulation. Environ. Model Softw. 35, 122–131. Calvari, S., Pinkerton, H., 1998. Formation of lava tubes and extensive flow field during the 1991–1993 eruption of Mount Etna. J. Geophys. Res. Solid Earth 103 (B11), 27291–27301. Camp, V.E., Roobol, M.J., 1989. The Arabian continental alkali basalt province: part I. Evolution of Harrat Rahat, Kingdom of Saudi Arabia. Geol. Soc. Am. Bull. 101 (1), 71–95. Camp, V.E. and Roobol, M.J., 1991. Geologic map of the Cenozoic lava field of Harrat Rahat, Kingdom of Saudi Arabia, . Saudi Arabian Directorate General of Mineral Resources Geosciences Map GM-123, scale 1:250,000: 1-37. Camp, V., Hooper, P., Roobol, M.J., White, D.L., 1987. The Madinah eruption, Saudi Arabia: magma mixing and simultaneous extrusion of three basaltic chemical types. Bull. Volcanol. 49 (2), 489–508. Cappello, A., Hérault, A., Bilotta, G., Ganci, G., Del Negro, C., 2015a. MAGFLOW: a physicsbased model for the dynamics of lava-flow emplacement. Geol. Soc. Lond., Spec. Publ. 426.

29

Cappello, A., Zanon, V., Del Negro, C., Ferreira, T.J.L., Queiroz, M.G.P.S., 2015b. Exploring lava-flow hazards at Pico Island, Azores Archipelago (Portugal). Terra Nova 27 (2), 156–161. Cashman, K.V., Thornber, C., Kauahikaua, J.P., 1999. Cooling and crystallization of lava in open channels, and the transition of Pāhoehoe Lava to 'A'ā. Bull. Volcanol. 61 (5), 306–323. Cashman, K.V., Soule, S.A., Mackey, B.H., Deligne, N.I., Deardorff, N.D., Dietterich, H.R., 2013. How lava flows: new insights from applications of lidar technologies to lava flow studies. Geosphere 9 (3), 1664–1680. Castruccio, A., Rust, A.C., Sparks, R.S.J., 2013. Evolution of crust- and core-dominated lava flows using scaling analysis. Bull. Volcanol. 75 (1), 1–15. Connor, L.J., Connor, C.B., Meliksetian, K., Savov, I., 2012. Probabilistic approach to modeling lava flow inundation: a lava flow hazard assessment for a nuclear facility in Armenia. J. Appl. Volcanol. 1, 3. Crisci, G.M., Rongo, R., Di Gregorio, S., Spataro, W., 2004. The simulation model SCIARA: the 1991 and 2001 lava flows at Mount Etna. J. Volcanol. Geotherm. Res. 132, 253–267. Deardorff, N., Cashman, K., 2012. Emplacement conditions of the c. 1,600-year BP Collier Cone lava flow, Oregon: a LiDAR investigation. Bull. Volcanol. 74 (9), 2051–2066. Del Negro, C., Fortuna, L., Herault, A., Vicari, A., 2008. Simulations of the 2004 lava flow at Etna volcano using the MAGFLOW cellular automata model. Bull. Volcanol. 70 (7), 805–812. Del Negro, C., Cappello, A., Neri, M., Bilotta, G., Herault, A., Ganci, G., 2013. Lava flow hazards at Mount Etna: constraints imposed by eruptive history and numerical simulations. Sci. Rep. 3 (3493). Dietterich, H.R., Cashman, K.V., 2014. Channel networks within lava flows: formation, evolution, and implications for flow behavior. J. Geophys. Res. Earth Surf. 119 (8) (2014JF003103). Dietterich, H.R., Soule, S.A., Cashman, K.V., Mackey, B.H., 2015. Lava Flows in 3D, Hawaiian Volcanoes. John Wiley & Sons, Inc., pp. 483–505. Egozi, R., Ashmore, P., 2008. Defining and measuring braiding intensity. Earth Surf. Process. Landf. 33 (14), 2121–2138. El Difrawy, M.A., Runge, M.G., Moufti, M.R., Cronin, S.J., Bebbington, M., 2013. A first hazard analysis of the Quaternary Harrat Al-Madinah volcanic field, Saudi Arabia. J. Volcanol. Geotherm. Res. 267, 39–46. El-Masry, N.N., Moufti, M.R.H., Nemeth, K., Murcia, H., Qaddah, A.A., Abdelwahed, M.F., 2013. Historical accounts of the AD 1256 eruption near Al-Madinah. VORISA Scientific Meeting, Jeddah, Kingdom of Saudi Arabia, pp. 9–13. Favalli, M., Pareschi, M.T., Neri, A., Isola, I., 2005. Forecasting lava flow paths by a stochastic approach. Geophys. Res. Lett. 32 (3). http://dx.doi.org/10.1029/2004GL021718. Favalli, M., Karátson, D., Mazzarini, F., Pareschi, M.T., Boschi, E., 2009. Morphometry of scoria cones located on a volcano flank: a case study from Mt. Etna (Italy), based on high-resolution LiDAR data. J. Volcanol. Geotherm. Res. 186 (3-4), 320–330. Favalli, M., Fornaciai, A., Mazzarini, F., Harris, A., Neri, M., Behncke, B., Pareschi, M.T., Tarquini, S., Boschi, E., 2010. Evolution of an active lava flow field using a multitemporal LIDAR acquisition. J. Geophys. Res. Solid Earth 115 (B11203). Felpeto, A., Araña, V., Ortiz, R., Astiz, M., García, A., 2001. Assessment and modelling of lava flow hazard on Lanzarote (Canary Islands). Nat. Hazards 23 (2-3), 247–257. Fries, C.J., 1953. Volumes and weights of pyroclastic material, lava, and water erupted by Paricutin volcano, Michoacan, Mexico. Am. Geophys. Union Trans. 34, 603–616. Ganci, G., Vicari, A., Cappello, A., Del Negro, C., 2012. An emergent strategy for volcano hazard assessment: from thermal satellite monitoring to lava flow modeling. Remote Sens. Environ. 119, 197–207. Griffiths, R.W., 2000. The dynamics of lava flows. Annu. Rev. Fluid Mech. 32, 477–518. Griffiths, R.W., Fink, J.H., 1993. Effects of surface cooling on the spreading of lava flows and domes. J. Fluid Mech. 252, 667–702. Guilbaud, M.-N., Self, S., Thordarson, T., Blake, S., 2005. Morphology, surface structures, and emplacement of lavas produced by Laki, A.D. 1783–1784. Geol. Soc. Am. Spec. Pap. 396, 81–102. Guilbaud, M.-N., Siebe, C., Layer, P., Salinas, S., Castro-Govea, R., Garduño-Monroy, V.H., Corvec, N.L., 2011. Geology, geochronology, and tectonic setting of the Jorullo Volcano region, Michoacán, México. J. Volcanol. Geotherm. Res. 201 (1-4), 97–112. Harris, A.J.L., Rowland, S.K., 2001. FLOWGO: a kinematic thermorheological model for lava flowing in a channel. Bull. Volcanol. 63, 20–44. Harris, A.J.L., Rowland, S.K., 2009. Effusion rate controls on lava flow length and the role of heat loss: a review. In: Hoskuldsson, A., Thordarson, T., Larsen, G., Self, S., Rowland, S. (Eds.), The Legacy of George P.L. WalkerSpecial Publications of IAVCEI 2. Geological Society of London, London, United Kingdom, pp. 33–51. Harris, A.J.L., Rowland, S.K., 2015. FLOWGO 2012, Hawaiian Volcanoes. John Wiley & Sons Inc., pp. 457–481. Harris, A.J.L., Murray, J.B., Aries, S.E., Davies, M.A., Flynn, L.P., Wooster, M.J., Wright, R., Rothery, D.A., 2000. Effusion rate trends at Etna and Krafla and their implications for eruptive mechanisms. J. Volcanol. Geotherm. Res. 102 (3-4), 237–270. Helz, R., Thornber, C., 1987. Geothermometry of Kilauea Iki lava lake, Hawaii. Bull. Volcanol. 49 (5), 651–668. Herault, A., Vicari, A., Ciraudo, A., Del Negro, C., 2009. Forecasting lava flow hazard during the 2006 Etna eruption: using the MAGFLOW cellular automata model. Comput. Geol. 35 (5), 1050–1060. Hidaka, M., Goto, A., Susumu, U., Fujita, E., 2005. VTFS project: development of the lava flow simulation code LavaSIM with a model for three-dimensional convection, spreading, and solidification. Geochem. Geophys. Geosyst. 6 (7), Q07008. Kauahikaua, J., Cashman, K., Clague, D., Champion, D., Hagstrum, J., 2002. Emplacement of the most recent lava flows on Hualālai Volcano, Hawai'i. Bull. Volcanol. 64 (3-4), 229–253. Kawabata, E., Cronin, S.J., Bebbington, M.S., Moufti, M.R.H., El-Masry, N., Wang, T., 2015. Identifying multiple eruption phases from a compound tephra blanket: an example of the AD1256 Al-Madinah eruption, Saudi Arabia. Bull. Volcanol. 77 (1), 1–13.

30

G. Kereszturi et al. / Journal of Volcanology and Geothermal Research 309 (2016) 14–30

Kereszturi, G., Németh, K., 2012. Structural and morphometric irregularities of eroded Pliocene scoria cones at the Bakony–Balaton Highland Volcanic Field, Hungary. Geomorphology 136 (1), 45–58. Kereszturi, G., Jordan, G., Németh, K., Dóniz-Páez, J.F., 2012a. Syn-eruptive morphometric variability of monogenetic scoria cones. Bull. Volcanol. 74 (9), 2171–2185. Kereszturi, G., Procter, J., Cronin, J.S., Németh, K., Bebbington, M., Lindsay, J., 2012b. LiDARbased quantification of lava flow susceptibility in the City of Auckland (New Zealand). Remote Sens. Environ. 125, 198–213. Kereszturi, G., Németh, K., Cronin, J.S., Agustin-Flores, J., Smith, I.E.M., Lindsay, J., 2013. A model for calculating eruptive volumes for monogenetic volcanoes — implication for the Quaternary Auckland Volcanic Field, New Zealand. J. Volcanol. Geotherm. Res. 266, 16–33. Kereszturi, G., Cappello, A., Ganci, G., Procter, J., Németh, K., Del Negro, C., Cronin, S.J., 2014. Numerical simulation of basaltic lava flows in the Auckland Volcanic Field, New Zealand — implication for volcanic hazard assessment. Bull. Volcanol. 76 (11), 1–17. Kerr, R.C., Lyman, A.W., 2007. Importance of surface crust strength during the flow of the 1988–1990 andesite lava of Lonquimay Volcano, Chile. J. Geophys. Res. Solid Earth 112 (B3) (n/a-n/a). Kerr, R.C., Griffiths, R.W., Cashman, K.V., 2006. Formation of channelized lava flows on an unconfined slope. J. Geophys. Res. Solid Earth 111 (B10), B10206. Keszthelyi, L.P., Pieri, D.C., 1993. Emplacement of the 75-km-long Carrizozo lava flow field, south-central New Mexico. J. Volcanol. Geotherm. Res. 59 (1-2), 59–75. Krauskopf, K.B., 1948. Lava movement at Paricutin volcano, Mexico. Geol. Soc. Am. Bull. 59 (12), 1267–1284. Macdonald, G.A., 1953. Pahoehoe, aa, and block lava. Am. J. Sci. 251 (3), 169–191. Moufti, M.R., Németh, K., 2013. The intra-continental Al Madinah Volcanic Field, Western Saudi Arabia: a proposal to establish Harrat Al Madinah as the first volcanic geopark in the Kingdom of Saudi Arabia. Geoheritage 5 (3), 185–206. Moufti, M.R., Moghazi, A.M., Ali, K.A., 2012. Geochemistry and Sr–Nd–Pb isotopic composition of the Harrat Al-Madinah Volcanic Field, Saudi Arabia. Gondwana Res. 21 (2–3), 670–689. Moufti, M.R., Moghazi, A.M., Ali, K.A., 2013. 40Ar/39Ar geochronology of the Neogene–Quaternary Harrat Al-Madinah intercontinental volcanic field, Saudi Arabia: implications for duration and migration of volcanic activity. J. Asian Earth Sci. 62, 253–268. Murcia, H., Németh, K., Moufti, M.R., Lindsay, J.M., El-Masry, N., Cronin, S.J., Qaddah, A., Smith, I.E.M., 2014. Late Holocene lava flow morphotypes of northern Harrat Rahat, Kingdom of Saudi Arabia: implications for the description of continental lava fields. J. Asian Earth Sci. 84, 131–145. Murcia, H., Németh, K., El-Masry, N.N., Lindsay, J.M., Moufti, M.R.H., Wameyo, P., Cronin, S.J., Smith, I.E.M., Kereszturi, G., 2015. The Al-Du'aythah volcanic cones, Al-Madinah City: implications for volcanic hazards in northern Harrat Rahat, Kingdom of Saudi Arabia. Bull. Volcanol. 77 (6), 1–19. Murcia, H., Lindsay, J.M., Smith, I.E.M., Moufti, R., Németh, K., Cronin, S.J., El-Masry, N.N., Niedermann, S. and McGee, L.E., under review. Petrogenesis of Late Quaternary basaltic magmas in northern Harrat Rahat, Kingdom of Saudi Arabia: Insights from whole-rock and mineral geochemistry, including noble gases in olivine Bulletin of Volcanology. Pellaton, C., 1981. Geologic map of the Al Madinah quadrangle. Ministry of Petroleum and Mineral resources. Geologic map GM-52C 250: Sheet 24D. Pinkerton, H., Wilson, L., 1994. Factors controlling the lengths of channel-fed lava flows. Bull. Volcanol. 56, 108–120.

Prewitt, J., 1970. Object enhancement and extraction. In: Lipkin, B.S., Rosenfeld, A. (Eds.), Picture Processing and Psychopictorics. Academic Press, Orlando, USA, pp. 75–149. Riker, J.M., Cashman, K.V., Kauahikaua, J.P., Montierth, C.M., 2009. The length of channelized lava flows: insight from the 1859 eruption of Mauna Loa Volcano, Hawai'i. J. Volcanol. Geotherm. Res. 183 (3–4), 139–156. Rossi, M.J., 1997. Morphology of the 1984 open-channel lava flow at Krafla volcano, northern Iceland. Geomorphology 20 (1-2), 95–112. Rowland, S., Walker, G.L., 1987. Toothpaste lava: characteristics and origin of a lava structural type transitional between pahoehoe and aa. Bull. Volcanol. 49 (4), 631–641. Rowland, S.K., Walker, G.P.L., 1990. Pahoehoe and aa in Hawaii: volumetric flow rate controls the lava structure. Bull. Volcanol. 52 (8), 615–628. Rowland, S., Garbeil, H., Harris, A., 2005. Lengths and hazards from channel-fed lava flows on Mauna Loa, Hawai'i, determined from thermal and downslope modeling with FLOWGO. Bull. Volcanol. 67, 634–647. Self, S., Keszthelyi, L., Thordarson, T., 1998. The importance of pahoehoe. Annu. Rev. Earth Planet. Sci. 26, 81–110. Soule, S.A., Cashman, K.V., Kauahikaua, J.P., 2004. Examining flow emplacement through the surface morphology of three rapidly emplaced, solidified lava flows, Kīlauea Volcano, Hawai'i. Bull. Volcanol. 66 (1), 1–14. Stoeser, D.B., Camp, V.E., 1985. Pan-African microplate accretion of the Arabian Shield. Geol. Soc. Am. Bull. 96 (7), 817–826. Swanson, D.A., 1973. Pahoehoe flows from the 1969–1971 Mauna Ulu eruption, Kilauea Volcano, Hawaii. Geol. Soc. Am. Bull. 84, 615–626. Tarquini, S., Favalli, M., Mazzarini, F., Isola, I., Fornaciai, A., 2012. Morphometric analysis of lava flow units: case study over LIDAR-derived topography at Mount Etna, Italy. J. Volcanol. Geotherm. Res. 235–236, 11–22. Thordarson, T., Self, S., 1993. The Laki (Skaftar Fires) and Grimsvotn eruptions in 1783–85. Bull. Volcanol. 55, 233–263. Valentine, G.A., Perry, F.V., 2006. Decreasing magmatic footprints of individual volcanoes in a waning basaltic field. Geophys. Res. Lett. 33 (14). http://dx.doi.org/10.1029/ 2006GL026743. Vicari, A., Alexis, H., Del Negro, C., Coltelli, M., Marsella, M., Proietti, C., 2007. Modeling of the 2001 lava flow at Etna volcano by a cellular automata approach. Environ. Model Softw. 22 (10), 1465–1471. Vicari, A., Ganci, G., Behncke, B., Cappello, A., Neri, M., Del Negro, C., 2011. Near-real-time forecasting of lava flow hazards during the 12–13 January 2011 Etna eruption. Geophys. Res. Lett. 38 (13). http://dx.doi.org/10.1029/2011gl047545. Walker, G.P.L., 1973. Lengths of lava flows. Philos. Trans. R. Soc. Lond. A Math. Phys. Sci. 274, 107–118. Wright, R., Garbeil, H., Harris, A.J.L., 2008. Using infrared satellite data to drive a thermorheological/stochastic lava flow emplacement model: a method for near real-time volcanic hazard assessment. Geophys. Res. Lett. 35 (19). http://dx.doi.org/10.1029/ 2008GL035228. Xu, W., Jónsson, S., 2014. The 2007–8 volcanic eruption on Jebel at Tair island (Red Sea) observed by satellite radar and optical images. Bull. Volcanol. 76 (2), 1–14. Zaidi, F.K., Mukhopadhyay, M., 2015. Morphometric analysis of the Scoria Cones and drainage pattern for the Quaternary and older volcanic fields in parts of the Large Igneous Province (LIP), Saudi Arabia. J. Afr. Earth Sci. 110, 1–13.