Volcanic Hazard Assessment of a Monogenetic Volcanic Field with Sporadic and Limited Information: Deterministic Approach for Harrat Lunayyir, Saudi Arabia

1. Introduction

Monogenetic volcanic fields are a common form of volcanism on Earth, especially within intracontinental regions [1,2,3]. They produce many small-volume volcanoes that usually erupt over short timescales, generating less than 1 km³ of dense rock equivalent (DRE) during their entire lifecycle [2]. Future eruptions typically occur at sites distant from previous vent locations. However, regions with ongoing melting anomalies may see later eruptions on top of older volcanic structures [4,5]. Fields with high magma output and frequent eruptions—ranging from centuries to thousands of years—can develop numerous volcanoes that overlap or form clusters, creating complex volcanic landscapes that might differ from typical monogenetic patterns [3].

When eruption recurrence intervals are geologically short, distinguishing individual eruptive units within such complexes becomes difficult due to the absence of erosional or sedimentological markers. These overlapping volcanoes may still be considered monogenetic if each is fed by a distinct conduit that solidifies before subsequent eruptions occur nearby. The spatial distribution of monogenetic vents varies widely, from random to clustered or structurally aligned patterns, often interpreted as reflecting subsurface melting anomalies and crustal structure [6,7].

Recent research utilizing machine learning suggests that mature volcanic fields can form massifs dominated by volcanic products, fostering vent clustering and potentially creating magma storage zones such as dykes and sills [8,9,10,11,12,13,14,15]. This can promote fractional crystallization and the formation of more silicic magmas, leading to unique explosive or dome-forming eruptions. Although these traits are commonly associated with polygenetic volcanism, in some instances within monogenetic fields, eruptions remain isolated events, forming solitary block-and-ash flows or lava domes. Evidence also indicates that monogenetic systems may gradually transition toward polygenetic behavior, with instances of multiple magma sources or magmatic evolution during brief eruptive sequences, complicating classification and study [16,17,18,19,20,21,22,23,24,25,26].

In summary, monogenetic volcanoes exhibit a complexity that exceeds their small size and brief eruption periods, reflecting their connections to various lithospheric processes occurring across different spatial and temporal scales. This complexity complicates the assessment of volcanic hazards in these areas. Key hazard considerations include the location and timing of potential future eruptions, as well as the nature of possible volcanic events. Effective hazard management requires probabilistic assessments based on data about the temporal and spatial patterns of volcanism, supported by comprehensive volcano models that can be incorporated into spatiotemporal prediction frameworks.

To build such models, three essential components are needed: a comprehensive event library with accurate geochronology spread evenly across the volcanic field; a detailed catalog of vents, including an understanding of vent distributions ideally placed within event chronology; and a simplified volcano model capable of capturing and synthesizing the sequence of volcanic activity observed at individual volcanoes. The accuracy of hazard estimates is thus directly affected by the quality and completeness of these parameters, requiring minimum threshold values for practical probabilistic approaches.

For example, if the event age database is incomplete, unevenly distributed, or characterized by significant uncertainties, whether due to methodological errors or challenges in establishing field-based stratigraphy, the reliability of probabilistic models is greatly diminished. Similarly, a limited understanding of vent distributions or the lack of a comprehensive vent catalogue can pose major obstacles to model development. An additional complication comes from the frequent occurrence of multi-event or nested volcanoes [27,28], where multiple eruptions form new edifices over relatively short timescales, further complicating the accuracy of probabilistic modeling in these systems.

A major challenge in developing probabilistic models is that eruption prediction timelines are often expected to be much shorter than what geological records support. This expectation is influenced by public familiarity with more frequent hazards, such as earthquakes, floods, or landslides. These events occur on timescales—ranging from years to centuries—that align with human planning horizons and are reflected in probability estimates for disaster preparedness. In contrast, volcanic eruptions—especially in association with distributed volcanic fields—typically have recurrence intervals far longer than 50–200 years, creating a paradox between expectations and reality [29,30,31,32,33,34,35,36,37,38,39,40,41,42,43]. No definitive solution exists yet. In this report, we offer a narrative to illustrate these challenges and suggest possible methods for addressing them elsewhere.

Saudi Arabia is experiencing rapid socio-economic transformation as it prepares for a sustainable, post-oil future. The country has many potential tourism sites, including recreation and adventure tourism, many of which require significant infrastructure improvements to accommodate visitors [44,45,46,47,48,49,50,51,52,53,54,55,56,57,58]. In recent years, the geoheritage of the region is more recognized [59,60,61,62,63,64,65,66]. However, additional, comprehensive geoheritage research is needed to identify key geological features and suitable geosites for inclusion in the national tourism strategy.

Significant development initiatives along the Red Sea coast focus on adventure tourism [67,68,69,70] amid the striking landscapes of the Arabian Shield, characterized by dark basaltic lava flows and relatively young volcanic formations sitting over light-colored Precambrian crystalline basement rocks. While Cenozoic volcanic fields are widespread across western Arabia (Figure 1 – harrats and Lunayyir Position – main roads and towns, etc.), current infrastructure limitations restrict the growth of niche tourism in the region. Harrat Lunayyir, situated approximately 300 kilometers northwest of Al Madinah and near Um Lujj, emerges as a key destination for specialized tourism activities (Figure 1). Red Sea Global plans to further develop tourism offerings in this area, leveraging the region’s dramatic contrast between basaltic lava fields and crystalline basement rocks [71,72].

During the planning stage, tourism development centered on the so-called Option 4 region (Figure 2 – overview of Lunayyir), which includes a young scoria cone (Target Volcano), its surrounding plain, and a complex lava flow that extends over 10 km from the cone. The region's general physiography is illustrated in Figure 3 and Figure 4, highlighting a significant morphological feature where Precambrian basement rocks emerge, creating a landscape of steep basement hills and flat-bottomed Quaternary valley networks following major structural elements [72]. These valleys are the geomorphological features that host most of the Cenozoic volcanics.

Geological information is compiled from the available 1:250,000 geological map of the region [GM-053 – Wadi al’Ays sheet] and an updated, interim dataset of ongoing 1:100,000 geological mapping. We differentiate lava flow horizons (and associated explosive pyroclastic eruptions) by identifying an older part of the area, labeled as the Qj lithostratigraphy unit, which is inferred to have formed in the Late Pliocene to Pleistocene and largely crops out in the southern part of the volcanic field. Younger volcanism from the Pleistocene to Holocene is distinguished based on sporadic radiometric age dating (with unreliable K-Ar and sporadic Ar-Ar data) [71,72,74,75,76,77,78,79,80,81,82,83,84,85], morphology [86], and geological mapping based on relative stratigraphy, labeled from Qm1 to Qm6 from older to younger (Figure 5 – 1 to 100 k geology map). The study area is primarily composed of the youngest Qm6 lava flows, with small remnants of Qm5 nearby and some older eruptive products that have been modified by erosion (Figure 5).

Preliminary research, however, shows that Option 4 is exposed to significant natural hazards, mainly seismic and volcanic risks [71]. Seismic data since 2007 (Figure 6 – volcanics and seismic), including a notable 2009 unrest [75,77,80,81,83,85], reveal ongoing low-intensity earthquake swarms along a NW-SE belt, related to shallow crustal structures and possible faults. The 2009 event—possibly a failed eruption—is valuable for anticipating future volcano-seismic activity and understanding warning signs of eruptions. Probabilistic and deterministic assessments indicate generally low earthquake intensity in the area, but similar events remain possible. The region also faces highly irregular hazards, such as rockfalls and flash floods. Option 4 lies in a zone with high volcanic vent density and includes one of the youngest volcanoes based on its youthful morphological appearance, suggesting potential for future eruptions even on a century scale.

The Option 4 area (Figure 5) is situated in the center of Harrat Lunayyir and has been designated as high-risk for both seismic and volcanic hazards, as well as potential cascading events such as mass movements. This understanding of the hazard profile at the proposed development site prompted an assessment that specifically concentrated on the volcanic risks in this region, especially within the Option 4 area. The evaluation was requested under two main constraints: 1) completing the report in approximately six months, and 2) providing suggestions regarding the likelihood of volcanic eruptions over a 50- to 200-year span, which aligns with planning needs for major investments and development projects.

While these requirements are generally achievable for hazards like earthquakes, floods, or landslides—thanks to established observation networks, historical records, and spatial analyses—they are more challenging for volcanic systems that have not been extensively studied. At Harrat Lunayyir, critical geological data are still unavailable; information on the age, composition, and eruption styles of the volcanoes is limited and unevenly distributed. This makes it difficult to compare with other Saudi Arabian volcanic fields, such as Harrat Al Madinah, which benefits from detailed geological maps and associated books [88]. Although interest in developing an eruptive history database grew after the 2009 volcano-seismic events, current data remains sparse and inconsistent, making probabilistic hazard estimates for the desired timeframes more complicated. To address these issues, a systematic approach was used to identify key volcanic hazards that require further attention. This report outlines the initial conceptual framework for the volcanic geology of the region, which was later developed into a detailed geological model. It also includes recommendations for risk mitigation and monitoring strategies as necessary.

2. Materials and Methods

This study uses extensive desktop research to identify global analogs of volcanism relevant to the study area, providing a set of comparative examples to describe potential volcanic hazards in the region. It incorporates previously documented field observations, new reconnaissance mapping, and field recordings along with macroscopic descriptions of key samples and microtextural observations from coherent (lava) and pyroclastic rocks (e.g., ash and lapilli). The methodology is divided into two main approaches.

The first step outlines the main types of volcanic hazards linked to the volcanic history of northern Harrat Lunayyir. These hazards are associated with specific eruption styles, as indicated by the volcanic deposits that can be mapped in the field.

Digital elevation models were used to calculate eruptive volumes, identify units, and characterize landscape features. The overall geomorphological assessment employed the ALOS-PALSAR Digital Elevation Model (DEM) with a 12.5-meter resolution, which is adequate for capturing general volcanic geological details and useful for lava flow inundation simulations with Q-LavHA software in QGIS environments [32,89,90,91]. Near the youngest vent (Option 4), LiDAR (Light Detection and Ranging) provided a processed DEM at a 0.5 m resolution from the Red Sea Global development company. Because LiDAR coverage was limited to the youngest and one older intact scoria cone, as well as the proximal part of the Holocene lava flow field, it does not allow for comprehensive high-resolution volume estimates across all eruptive products. Nonetheless, LiDAR data enabled a detailed analysis of lava flow surface textures, pressure ridges, lava tubes, and channels, and facilitated a comparison between volume estimates from available LiDAR field data (a small area only around the source volcano) and ALOS-PALSAR DEMs.

Red Relief Surface Maps (RRSM) [https://www.rrim.jp/en/ - accessed on 30 October 2025] were generated by combining Red Slope Angle Maps and Differential Openness maps derived from the DEMs [92,93]. RRSM proved useful for delineating lava flow boundaries, enabling the identification of individual lava lobes [24,94,95,96,97] with LiDAR data, while ALOS-PALSAR provided broader lava flow distribution patterns.

Eruptive and edifice volumes were calculated from the DEMs using the QGIS Volume tool, which operated with the lowest and average pixel elevations along identified lava flow contacts. Field observations and GPS elevation measurements revealed that most emplacement surfaces are remarkably flat, dipping less than 1.5 degrees (approximately a 300 m drop over 12,000 m). This slight dip is unlikely to significantly impact volume calculations, except in areas where topography drops sharply—mainly along suspected structural lines where lava flows cascade. These zones were identified in the distal parts of the study area, where flows narrow and thin, contributing minimally to the total estimated volumes.

3. Results

This section may include subheadings. It should give a clear and concise summary of the experimental results, their interpretation, and the conclusions that can be drawn.

3.1. Field Observations

3.1.1. Spatter and Scoria Cones

Fire or lava fountaining eruptions happen when jets of low-viscosity magma shoot molten fragments, called spatters, into the air [98,99]. These sometimes-episodic events, which can look like water hose blasts, last from hours to days and are usually linked to fissures hundreds of meters long before focusing into point-source eruptions. Fountain heights often go beyond 200 meters, so a 500-meter exclusion zone is advised to stay clear of hot pyroclasts [98,99]. Magma speeds exit at several tens of meters per second. Fissure eruptions produce lower fire curtains, with mass eruption rates ranging from 10⁵ to 10⁶ kg/s. Main products include large spatters around vents, forming cones, and new lava flows that can cover older volcanic features. In the lava fields west of Target Volcano, volcanic spattering has covered previous landforms (Figure 7).

Weak convective plumes above magma fountains only lift fine ash particles a few hundred meters, while larger fragments quickly fall out [98,99]. Fountaining often occurs alongside or shifts into effusive activity, creating interconnected lava flows and mild explosive hazards. Lava fountains typically feed flows both directly and by depositing spatter clasts, resulting in alternating agglutinates and clastogenic lava where clastic textures remain visible. These features are standard near the source area. Hawaiian fountaining involves evenly bubbly magma ascent, contrasting with Strombolian eruptions, which feature slug flow, large bubble bursts, and greater pyroclast dispersion [99].

Volcanoes that mainly emit mafic, low-silica magmas are characterized by periods of open vent activity. During these times, explosive eruptions typically occur in rapid, individual bursts of spatter—called Strombolian-style (sensu stricto) explosive events—separated by intervals ranging from seconds to several hours. Explosion durations typically last from a few seconds to several tens of seconds [99,100,101]. These events produce rapid sprays of mostly fluidal spatter, scoria, volcanic bombs, and small amounts of ash. The deposits build up within a few hundred meters of the vent, gradually creating landforms known as scoria (cinder) cones [102].

Due to localized pyroclast dispersion during typical Strombolian eruptions, most eruptive products directly contribute to cone development. From a volcanic hazard standpoint, Strombolian-style explosions can disperse pyroclasts up to 1000 meters from the eruption site; this distance provides a helpful guideline for exclusion zones near active scoria cones. These eruptions mainly produce spatter cones, ramparts, and structures made of coarse scoriaceous lapilli and bomb-rich deposits. When smaller but sustained eruption columns occur, larger scoria cones composed of fine lapilli and coarse ash can form, with localized fallout deposits building up within approximately one kilometer. Strombolian eruption columns or spray heights typically range from a few tens of meters to under 200 meters, ejecting several hundred cubic meters of pyroclastic material per burst. Pyroclast exit velocities vary widely from 3 to 400 m s−1 [99]. Spatter bombs and scoria lapilli are projected along ballistic, mostly parabolic trajectories, contributing to the formation of modest spatter and scoria cones or ramparts beside the vent. Minor amounts of ash can get caught in short-lived, weak convective plumes, sometimes carried away by prevailing winds, although such plumes rarely reach above a few kilometers [99].

The study area exhibits numerous examples of low-intensity Strombolian-style explosive events, which are often part of larger eruption sequences characterized by occasional higher-intensity explosions and the gradual development of scoria cones. The cone-building pyroclastic succession at the Target Volcano in Option 4 illustrates these deposit types.

Figure 8. The Target Volcano (25°20'33.69"N, 37°41'5.89"E), with its young scoria cone morphology (Qm6 vent), is visible in this view due to Strombolian-style explosive eruptions. An older scoria cone (25°22'21.57"N, 37°40'5.94"E), located at the left front of the Target Volcano, is part of the oldest volcanism in the volcanic field (Qm1). Its preserved vent is not as steep as the Target Volcano's vent.

Figure 8. The Target Volcano (25°20'33.69"N, 37°41'5.89"E), with its young scoria cone morphology (Qm6 vent), is visible in this view due to Strombolian-style explosive eruptions. An older scoria cone (25°22'21.57"N, 37°40'5.94"E), located at the left front of the Target Volcano, is part of the oldest volcanism in the volcanic field (Qm1). Its preserved vent is not as steep as the Target Volcano's vent.

3.1.2. Lava Flows

The formation of a solid crust on lava flows depends on how cooling influences viscosity over time, compared to the lava's shear rate [103]. In low-viscosity basalts, such as pāhoehoe, the viscoelastic layer absorbs strain, allowing the smooth crust to remain intact even when lava beneath flows rapidly, sometimes exceeding 10 m/s. This crust also remains preserved at low strain rates, such as when lava moves slowly across flat valley floors west of the Option 4 area. When strain rates on developing crust increase—due to steep slopes, sudden changes in slope angle, higher lava input, or increased eruption rates—the crust experiences plastic failure, resulting in surface textures that signal complex brittle-elastic deformation. This process leads to clinker, slabby, and ultimately aa-type lava flows (Figure 9).

As lava viscosity increases in distal regions and flow slows, its crystallinity rises, causing it to resemble ‘a‘ā lava [103]. Transitional surface textures often appear between smooth pāhoehoe and clinker-rich a’a, creating hybrid forms like shelly, slabby, and rubbly pāhoehoe—commonly seen in the study area (Figure 10) [104,105]. These transitions are also typical of the harrats in western Saudi Arabia [106].

When lava crusts form over a channel and stay intact, they create an underground lava tube, offering an alternative to open lava channels. In the Option 4 area, collapsed tube roofs reveal entry cavities (Figure 11), and nearby depressions, up to 10-15 meters deep, indicate extensive tube networks along the main systems. Lava tubes are also observed in other parts of Harrat Lunayyir, featuring short, well-developed tubes with visible openings in the central region, as well as near vents within the Option 4 area. Usually, tubes and open channels form complex distributary systems, found both near and away from volcano vents, especially where terrain slopes change. Smooth lava surfaces with fractured plates, as seen on satellite images, suggest the presence of complex tube networks extending several meters underground.

Additionally, two main open channels are identified. As cooling increased viscosity, the lava transformed into substantial rubble piles. Mechanical abrasion along these channels broke larger pāhoehoe crust plates into smaller, semi-rounded boulders. The marginal areas, which experienced more cooling and higher viscosity, formed toothpaste-like flows and characteristic 'aa' type lavas with squeeze-up features. The rough-surfaced region between the main open lava channels primarily consists of rubbly pāhoehoe. It is believed to be a former ponded lava zone, where initial crust formation was interrupted by new influxes of fresh lava, resulting in a complex rubbly pāhoehoe field. When a viscoelastic crust thickens along a lava channel, it solidifies, forming a strong roof over the flow. Molten lava then moves beneath this crust, keeping heat and mobility for longer periods. Large lava tubes, such as those in Harrat Rahat, can extend over a kilometer and serve as underground passages, allowing lava to travel great distances. The roof may collapse if the lava supply decreases or after eruptions, when lava drains away, leaving partially empty tubes. Fluctuations in lava output can add to surface complexity, blending collapsed tube roofs with active molten channels and viscous flows. This results in rugged terrain, as seen at Option 4’s proximal lava fields. Skylights over lava tubes form during eruptions, offering views of moving molten lava; after cooling, they become deep cavities linked to tube networks. These tubes, ranging from less than 1 m to over 30 m in diameter and 15 m in height, can channel lava over hundreds of kilometers beneath a mostly solid surface [107]. Surface features, such as pāhoehoe slabs, skylights, and updoming, indicate meter- to tens-of-meters-wide tubes along major lava pathways. Large roof collapses and extended tube openings suggest long, narrow tubes formed where lava accumulated behind barriers or in depressions. Lava stored in these main zones later traveled farther than expected. Tube-fed pāhoehoe flows can extend much farther than comparable a’a flows, making tube formation a key consideration in assessing lava flow hazards. A solid, insulating crust in pāhoehoe lavas enables flow inflation as magma accumulates beneath the rigid surface when the volumetric flow rate exceeds the advance rate, often at slope breaks [108,109]. This results in surface uplifts and can cause lava breakouts. Although most common in pāhoehoe, inflation and breakouts also occur in a’a lavas when their fronts slow due to cooling or gentler slopes. The resulting autobreccia increases resistance and thickens the flow. Pressure buildup behind stationary fronts may trigger bright breakouts through the margins or surface, a scenario typical in the study area and a hazard to consider during eruptions and recovery efforts.

3.1.3. Ash Plains

Explosive eruptions from low-viscosity magma, as observed at Mexico’s Parícutin volcano (1943–44), can produce large amounts of fine ash [110]. The initial phase of Parícutin's eruption began on February 20, 1943, quickly forming a 167 m scoria cone with up to 18 explosions per minute and a plume roughly 5 km high. Parícutin serves as a valuable comparison for the Target Volcano in Option 4 because of its similar sequence and structure. Between March 18 and June 9, 1943, the eruption intensified, creating more ash and a convective plume reaching 6 km and carried downwind [110]. The eruption style was characterized by frequent explosions, with a magma eruption rate estimated at 10⁴–10⁵ kg/s. Ash deposits spread as far as 250 km from the volcano, resulting in continuous tephra coverage; lapilli-sized fragments were transported several kilometers from the vent [111]. Inevitable scoria fall and cone-building deposits with high fragmentation and dispersal rates are commonly defined as "violent strombolian" [112,113], highlighting the episodic nature of these explosions despite differences from classic strombolian eruptions. Extensive ash and lapilli fields within Option 4, associated with the Target Volcano, suggest that violent Strombolian explosive eruptions have occurred in the region, indicating a potential volcanic hazard.

Recent assessments suggest using the term "violent Strombolian" specifically for eruptions that produce columns typically under 8 km and generate sizable fallout deposits with limited stratification [110]. However, vertical grain size variations can occur, as seen in ash and fine lapilli deposits approximately 3.5 km from the Target Volcano (Figure 12). These eruptions are characterized not only by their column height but also by sustained, pulsatory explosive activity that creates cone-shaped, stratified deposits rich in fine ash. In contrast, micro-Plinian eruptions tend to be more continuous and produce larger, more extensive deposits, such as coarse scoria cones and widespread massive deposits [98]. Recently, some eruption styles previously categorized as Strombolian have been reclassified as micro-Plinian, based on their consistent ash and pyroclast output [98]. This distinction highlights differences from violent Strombolian scoria cones. To better understand these styles, documented examples, like a detailed study of a scoria cone-forming eruption, are useful analogies. For instance, the 1973 basaltic eruption at Heimaey, Iceland, started with lava fountains reaching 50–100 meters from multiple vents along a 1.5 km fissure. [114,115]. Activity then focused on three vents near the fissure’s northern end, forming a 200-meter-high scoria cone with associated lava flows. Individual explosions resembled discrete fire fountain sprays. Later phases involved more sustained activity, producing a buoyant eruption column and depositing coarse, massive scoria near the vent, as well as localized fall deposits downwind. This event provides an analog to interpret potential eruption types and scales for the Target Volcano, informing future scenarios in the study area.

3.2. Volcanic Hazard Types

After outlining major eruption styles, we identify key volcanic hazards for the area. Using analogies and observations, we summarize these hazards and suggest mitigation measures in a table. We describe the precursory phase as the time of geophysical warnings, followed by hazard-specific onset stages, established eruptive processes, and post-event outcomes.

Large volumes and widespread lava flows indicate the importance of volcanic hazards associated with lava effusion. These can be classified by lava characteristics, such as pāhoehoe or transitional flows evolving into aa-type lava. Such hazards usually move slowly and are considered manageable, but they can cause significant long-term landscape changes and result in destruction.

Volcanic hazards can also arise from a range of explosive eruption styles, including: 1) the formation of lava spatter cones with associated lava fountaining; 2) the development of scoria cones; 3) Strombolian or mafic sub-Plinian explosive eruptions; and 4) deposition of ballistic projectiles. Additionally, volcanic degassing and vog (volcanic gas fog) are often observed alongside extensive lava flow emplacement, especially after the formation of scoria cones. Based on field observations and considerations, 8 different types of volcanic hazards can be identified. Each of them is documented in Appendix 1 in detail, such as

1)

Pāhoehoe lava flow outpour

2)

Transitional lava flow outpouring

3)

Lava spatter cone formation, lava fountaining

4)

Scoria (cinder) cone birth and growth

5)

Ballistic impact

6)

Cone collapse, cone rafting

7)

Sub-Plinian explosive eruption and associated ash fall

8)

Volcanic degassing, vog formation

3.3. Eruption Scenario Simulations

Eruption scenarios are concise, realistic narratives that summarize potential eruption types, timelines, and impacts, facilitating clear information transfer for volcanology, especially when probabilistic methods are not practical. Due to limited geological data at Harrat Lunayyir, hazard-based models were created using direct observations and desktop analysis. Two scenario types—localized point eruptions and fissure eruptions—have both occurred over the past 600,000 years, with fissure eruptions being more common and likely in future events. Three potential eruption sites were identified: vent 1 NW outside Option 4 (Figure 13); vent 2 near Target Volcano (Figure 14); and vent 4 SW inside Option 4 (Figure 15). Fissure vents were located at the same sites previously considered for point source vents (Figure 16 & Figure 17), with additional fissure settings simulating a particularly long fissure in the southwest of the region (Figure 18), aligned with the main structural feature of the area.

3.3.1. Well-Localized Eruption Point Development

Eruption points, typically singular locations, often form volcanic features like spatter or scoria cones. In our study, we used a 500 m radius for eruption footprints and a 1000 m exclusion zone for ballistic impacts. Lava flow simulations used the Q-LavHA QGIS plugin, with input parameters including a conservative 15 m lava thickness and a 25 m buffer, which provides a minimum estimate of the affected area. Any location with a non-zero probability in the simulation is considered at hazard. Because the terrain is flat, a flow of 50 km was set, following Q-LavHA recommendations for flat surfaces [89,90], consistent with known flows in Harrat Lunayyir. The Manhattan technique was applied to follow directional valleys typical of the region. Results are probability maps (Figure 13, Figure 14 & Figure 15) that mark all areas reached by lava as hazardous, with specific probabilities not emphasized.

3.3.2. Localized Fissure Eruption

After simulating eruptions from point sources, further tests examined how fissure eruptions could alter scenarios and impact regions (Figure 16,Figure 17&Figure 18). This was important for hazard assessments because most volcanoes in the Lunayyir field formed through fissures of different lengths. Short fissures, approximately a hundred meters long, often form elongated structures, while longer fissures extend several kilometers and feature aligned vents. One test modeled a 15 km NW-SW-trending fissure, like the 2009 event, with lava effusion simulated by Q-LavHA using emission points spaced every 90 m along its length (Figure 18). Each model included a 500 m buffer to represent areas where lava spatter cones and fountaining could deposit material near the active vents.

3.4. Lava Flow Inundation Hazard Zonation

This summary highlights the main lava inundation models (Figure 19), excluding cone base and ballistic impact visualizations, as their effects are much less significant. We used the lava inundation zone as a proxy to identify areas where development would be at high risk, which also aligns with previously recognized vent clusters, alignment zones, and graben features marked by structural weaknesses.

3.5. Lava Flow Volume Estimate

Estimating lava flow volume is key to understanding volcanic hazards. However, in our study area, there were no features suitable for accurate stratigraphic assignment or age dating due to time and resource constraints. As a result, we relied on existing maps, limited remote sensing, and targeted field surveys to outline the most probable extent of a young lava flow in northern Harrat Luinayyir. Comprehensive mapping of the nearly 100 km² region was not feasible within the project timeline, so surveys focused on key locations. Our findings, along with previous lava flow maps, showed varying interpretations of the youngest lava effusion in the area (Figure 20 & Figure 21).

The flow boundary delineation discrepancies were more apparent in the distal regions, where the distinction between main flow arteries, channels, and thin squeezed-out lava surfaces is difficult to distinguish. Red Image Surface Maps and a DEM combination map have also been used based on ALOS-PALSAR DEMs. As the DEM resolution is only 12.5 m, the resulting maps were only helpful in seeing the internal textures of the main lava bodies (Figure 22 & Figure 23), confirming the separation of lava sections within the flow field and allowing for slightly more restricted lava flow coverage, as previous mappings indicated.

For the lava volume calculation, the emplacement surface was generated using elevation data from around the flow margins and linear interpolation in QGIS. This surface is generally flat, with step-like drops in elevated areas and a gentle westward slope (Figure 24). The estimated thickness map indicates that the flow margins are thin, and the interior rarely exceeds a few tens of meters in thickness. The relative thickness map highlights key accumulation zones within the field (Figure 25). However, due to DEM resolution limits and uncertainties in pre-eruptive surface location, these values are approximate and used mainly to demonstrate that the flow followed the microtopography and ponded east of the cone (Figure 25). Assuming a uniform thickness of 10 m and an area of 30.616 km², the total lava volume is approximately 0.3 km³. Using an average lava thickness of about 30 m and uniform coverage, the total erupted lava volume would be around 1 km³. When accounting for voids, lava tubes, and high vesicularity, the dense rock equivalent (DRE) is approximately 30% of this volume, resulting in a range of 0.09 to 0.3 km³, based on estimates of 10 or 30 m lava thickness, respectively. These values are comparable to those of the well-documented Paricutin (Mexico) eruption [116], where the average flow thickness estimates ranged from 8 to 25 meters, while the covered area was calculated to be in the range of 1.892 to 12.317 km², yielding erupted volumes of 0.01 to 0.3 km³ without DRE corrections. DRE corrections would lower the eruptive volumes by approximately 20–30%. As Paricutin's activity occurred between February 1943 and March 1952 in 23 distinct eruptive phases, we can assume that the formation of the Qm6 lava flow field likely took place in multiple stages and could have lasted several months or more. It appears that our study area has not exhibited such complex and numerous eruptive episodes (23) and is likely to have experienced only a few, based on the lava flow mapping so far. This suggests that the eruption rate there might be estimated as higher than that of Paricutin. In short, in a shorter period with fewer eruptive phases, a larger volume of lava was produced. Volume calculations using GIS-based volume counts appear to provide overestimates of values ranging between approximately 1.6 and 3.3 km³, reflecting the uncertainty in defining the emplacement surface, the GIS-based calculation methods, and the resolution (12.5m) of the available DEM.

3.6. Ash Fall Hazards

Ash fall hazards in Saudi Arabia's Quaternary volcanic fields receive less attention than lava flow risks, mainly because ash deposits—often fine lapilli to fine ash—are rarely well-mapped. This is due to their preservation as mixed aeolian and pyroclastic layers, making mapping challenging. However, recent studies, especially near Al Madinah, have improved the understanding of these ash plains. The 1256 CE eruption southeast of Al Madinah produced a chain of scoria cones and a significant ash plain [117], similar to that found within the Option 4 area (Figure 26).

Based on limited spot measurements (due to access or preservation potential), a rough isopach map of tephra distribution has been created (Figure 27).

The United States Geological Survey (USGS) uses Ash3D software to predict the spread (ash clouds) and deposition (ash fall) from volcanic eruptions, primarily during Plinian eruption conditions (USGS Ash3D). This three-dimensional model simulates volcanic ash movement by factoring in time-varying wind patterns and other weather data to estimate ash transport under both current and past atmospheric conditions. The online tool provides near-real-time forecasts of ashfall centered on Harrat Lunayyir, using data from the Smithsonian Global Volcanism Program Holocene Volcanoes catalog. The vent is located near the center of the volcanic field, close to recent eruption sites southeast of the Option 4 area, roughly 20 km away. The model considers the estimated magma volume involved—specifically the tephra—and the height of the eruption column above sea level. Harrat Lunayyir sits about 1300 meters above sea level, so an eruption column reaching roughly 8700 meters above sea level corresponds to a 10 km plume in the model. This matches sub-Plinian or violent Strombolian eruptions. An initial simulation was run for an ash-producing eruption on July 19, 2024, using real-time wind data, with a plume height of 7 km above sea level (around 5.7 km above the volcano) and a total eruptive volume of 0.005 km³ DRE, with 5% of the volume airborne at the vent (Figure 28). The results, shown in Figure 31, suggest that central Harrat Lunayyir would receive between 1 and 3 centimeters of ash over an area three times larger than in Option 4. These modeled values are much lower than observed in the field and might reflect only a single phase of the eruption rather than the total effect across all eruption stages.

Ash fall hazard was identified during field mapping, with primary ash layers approximately 2 mm thick located about 8 km northwest of the youngest cone (Target Volcano). This distribution pattern indicates that a Strombolian or sub-Plinian style explosive eruption dispersed ash through the atmosphere over this distance. Based on these findings, an average ash distribution pattern can be constructed, assuming eruptions occur along the central axis or along a fissure aligned with the observed structural trend in the field. This approach provides a visual representation of the maximum extent, whereas the ash fall hazard should be considered. Since the impact of this hazard is generally low, it is recommended to acknowledge its presence, especially in spatially restricted areas where its effect could be more significant (such as potential roof collapses). Figure 29 illustrates this zone.

Ash fall can significantly impact buildings and their critical operational components. Inside, ash contamination may pose health risks to occupants (www.ivhhn.org), directly damage sensitive equipment such as electronics, and cause abrasion damage to flooring surfaces. Abrasive damage to roofing, cladding, and flooring materials during ash removal is often unavoidable, even with minimal ash buildup, which increases recovery costs. Ash fall may also disrupt essential building services, including electricity and water systems. Construction planning should consider using appropriate non-structural elements to reduce ash fall risks. Gutter systems, when installed, can collect ash from roofs, reducing drainage capacity and adding to structural loads; blockages in gutters and downpipes can lead to localized flooding and damage, especially to roofs, drainage networks, and ceiling spaces. Internal gutters are particularly vulnerable and often hard to clean. Metal roofs, fastenings, and claddings may be prone to corrosion, as ash deposits containing salt can accelerate deterioration. Painted surfaces and artistic features are at risk from acid leachates, which can cause rapid and significant damage, especially for tourism facilities that depend on aesthetic appeal. Structural concerns include potential roof deformation and collapse due to excessive ash loading. However, very thick deposits (>100 mm, more often >300 mm) are uncommon except within 5 km of the vent, they remain a design consideration. Field measurements recorded ash thicknesses of 300 mm at 4 km and 10 mm at 8 km from the source. Long-span, low-pitched roofs should generally be avoided, as they are most vulnerable to ash loading. Assuming an average basalt density of 2,900 kg/m³ and approximately 75% vesicularity observed in the study area, roof designs should be able to withstand up to 500 mm of ash cover (about 362 kg/m²), which is relevant within 2 km of a vent where lava inundation could occur within days. While simultaneous rainfall and ashfall are unlikely in this arid environment, the possibility cannot be ruled out. The overall impact level depends on these interconnected factors, as buildings function as integrated systems. Mitigation strategies, such as prompt and regular ash removal when safe, can lessen or even prevent adverse effects. Cleaning building exteriors and interiors after ashfall can be costly and time-consuming and may require ongoing effort due to recurring ash or wind- and water-driven remobilization. Taking key steps before, during, and after ashfall can help minimize damage. Protective gear is essential during cleanup, as roofs may be slippery or collapsing, especially if water blasting increases the ash’s weight. Fine ash removal creates dust that can cause discomfort, and all collected ash should be disposed of where it won’t be remobilized later.

4. Discussion

Geological studies and simulations indicate that Option 4 is in an area highly vulnerable to volcanic hazards, particularly lava flows, with eruptions potentially reaching a VEI 4 severity and dispersing ash over tens of kilometers. The likely eruption patterns closely resemble those seen in Iceland over the last thousand years, due to similar rift-related processes, basaltic sources, and fissure vent formations. Recent well-documented eruptions in Iceland have produced spatter cones, mild explosive phases, and steady low-viscosity lava flows, offering valuable insights for hazard mitigation in Harrat Lunayyir.

Volcanic hazards in the Option 4 area and surrounding regions are classified as having a low to moderate risk, although they can still cause significant impacts. Finding safe zones within this region is challenging, as approximately 80% of the Option 4 area is exposed to one or more types of volcanic hazards (Figure 30). In addition to lava flows, ash fall poses a significant threat, and an eruption could potentially reach VEI 3 or VEI 4 levels, producing extensive ash plains that may block access and mobility. Ash fall also increases the risk of total loss of facilities, roads, and equipment planned for development and use in Option 4.

There is currently no established industry practice for mitigating volcanic hazards related to monogenetic volcanism, especially in Saudi Arabia, where modern volcanic eruptions have not been observed. The most recent confirmed eruption occurred in 1256 CE southeast of Al Madinah City [118]. Stratigraphic evidence also suggests that eruptions within Harrat Lunayyir may date back to the Holocene, with an estimated occurrence between 4,500 and 1,000 years ago at the Target Volcano, which could be among the youngest sites in the region. However, these age estimates are uncertain, and the likely recurrence interval probably falls well outside a typical 50-year timeframe. This estimate is based on the number of vents—usually five to ten per volcanic complex—leading to an average recurrence rate of decades to centuries. However, this method highlights the limitations of current models.

Volcanic hazards in monogenetic volcanic fields—characterized by small volumes, mainly basaltic composition, absence of central magma chambers, dispersed volcanoes over large areas, vent distribution dictated by geology, and brief eruption durations—while the overall activity can last hundreds of thousands of years—pose significant challenges from a probabilistic standpoint. Often, there is limited age data or large error margins, making it difficult to predict where and when volcanoes form and erupt accurately; these two aspects are distinct, as shown by Harrat Lunayyir, where many volcanoes are complex vent structures that can become active rapidly but are supplied by different magma sources from separate melting events.

Industry practices have generally addressed volcanic hazards from monogenetic volcanism only at high-sensitivity sites, such as nuclear power plants, dams, major infrastructure, or high-risk waste repositories, where design criteria must consider hazards identified over thousands of years. For instance, high-level waste repositories (usually for nuclear waste) require hazard assessments spanning about 50,000 years—a period during which multiple vents can form in volcanic fields with eruption recurrence rates over thousands of years, as seen in the Auckland Volcanic Field (New Zealand) [37,119,120,121], Yucca Mountain (Nevada) [12,122,123], and Southern Armenia [35,124,125]. Studies in these regions have focused on urbanization and lifeline exposure (Auckland), planning for high-level power plants (Armenia), or designing waste repositories (Yucca Mountain). These studies were conducted within appropriate time frames and supported by high-resolution, evenly distributed age dating across the fields.

In contrast, predicting volcanic hazards for Harrat Lunayyir over a 50-year period is currently not feasible due to data and technical limitations. Even the most precise Ar-Ar radiometric dating methods have errors that surpass the resolution needed for this timescale. Consequently, the approach shifted from probabilistic estimates to a deterministic method that systematically considers all identified volcanic hazards. This involves preparing for any documented hazards throughout the entire field (~600,000 years of records), with a particular focus on Qm5 and Qm6 vent locations. These sites span an area like the broader Lunayyir vent zone, supporting the strategy of evaluating hazards based on the overall vent distribution of the field.

For industry applications, each identified type of volcanic hazard should be considered a potential scenario at any location where vents have historically formed. Design solutions should be scaled accordingly: for example, structures need to withstand lava flow inundation by considering fire hazards, mechanical stress, and assessing potential economic losses from such incidents; roofs should be designed to support at least 10 cm of ash accumulation, ideally up to 50 cm; and precautions should be implemented to protect sensitive equipment from fine ash infiltration, such as air conditioners, vehicle and generator filters, sewage systems, and drinking water supplies.

Hazard zones were identified based on three simulated long fissure lava flows, using the orientation and length of the 2009 fissure (~15 km) as a model. Using Q-LavHA, vents were placed every 100 meters along each fissure, forming three NW-SE trending lines (west, central, east). The lava inundation maps were combined and smoothed with a 100-meter buffer to define the ORANGE zone, indicating moderate volcanic hazard from lava flows. The RED zone, representing the highest hazard, was established by adding a 1,000-meter buffer around the fissures, where immediate edifice formation and pyroclastic activity are primary threats. YELLOW zones denote lower risk areas mainly affected by ash fall up to 1 cm, as shown in Figure 25.

Figure 31. A simplified hazard zonation map derived from lava flow inundation simulations based on three extended fissure scenarios. The fissure lengths reflect the parameters observed during the failed 2009 eruption and subsequent fissure activity. Zone 2, or the medium hazard zone, represents the combined inundated areas from all three simulations. The primary NW-SE structural trends in the region, documented vent alignments, and the broader geological framework of greater Harrat Lunayyir informed the selection of fissure locations.

Figure 31. A simplified hazard zonation map derived from lava flow inundation simulations based on three extended fissure scenarios. The fissure lengths reflect the parameters observed during the failed 2009 eruption and subsequent fissure activity. Zone 2, or the medium hazard zone, represents the combined inundated areas from all three simulations. The primary NW-SE structural trends in the region, documented vent alignments, and the broader geological framework of greater Harrat Lunayyir informed the selection of fissure locations.

High-hazard zones are defined as the area immediately surrounding each fissure, including a 2,000 m buffer to account for potential proximal hazards such as ash fall, ballistic projectiles, lava fountain reach, and scoria cone growth. Low-hazard zones extend up to 8,000 m from the identified fissures, where an ash accumulation of at least 1 cm is anticipated. This threshold is supported by field studies and mapping, serving as a reliable indicator of potential volcanic ash hazard, especially if a future sub-Plinian eruption occurs—a scenario with significant likelihood based on the area's geological history.

It should be noted that the southeast corner of the study area corresponds to the site of the 2009 fissure opening and features a historically high vent density and explosive sub-Plinian activity. For the hazard map, volcanic hazards related to this southeastern region have not been included, allowing for a focused assessment of the northern territory. The north area, where Option 4 is located, remains clearly separated from central regions by prominent basement horsts, which influence fissure propagation and potential vent formation. These morphological features also act as practical barriers to prevent lava flow inundation from the southeast into the Option 4 area. Although ash fall may still pose a significant hazard, it does not materially change the overall hazard classification for Option 4. Additionally, this map does not include information on secondary or cascading hazards, such as remobilized ash or potential flash flood routes, nor does it address extremely distal, thin lava flow inundation. Thin lava flows are considered low hazard, while comprehensive modeling of flash flood risks requires further research focused explicitly on sedimentary gravity flows. Generally, flash flood channels or alluvial fans are usually located along steep valley slopes or near the foothills of major mountain ranges.

Estimating eruption probability remains uncertain at this stage due to limited precise knowledge about volcanic event sequences, which mainly relies on approximate assessments based on relative stratigraphy. The total activity of the volcanic field spans approximately 600,000 years, encompassing over 700 vent openings across around 150 multi-vent volcanoes. This age estimate is supported by 17 high-precision Ar-Ar radiogenic age dates, mostly from older lava flows, though the dataset is geographically and temporally uneven [84]. Previously, the total lifespan of the volcanic field was estimated to be 2.8 million years, which adds further uncertainty to probability calculations. The average recurrence rate for vent openings over the 600,000-year period would be roughly one event every 857 years, which can be rounded to a 1,000-year recurrence interval if events were evenly spaced—although this uniform distribution is unlikely. Therefore, the 1,000-year vent-opening interval should be viewed as a conservative estimate, with actual recurrence rates likely being higher. When considering the formation of volcanic edifices, with an estimated 150 volcanoes, the recurrence rate would be about 4,000 years. Based on these figures, it is reasonable to expect that complex, multi-vent eruptions could happen roughly every 4,000 years, as indicated by geological mapping throughout the entire field, regardless of vent age or location.

Assuming Qm5 began 4500 years ago and Qm6 started 1000 years ago, there were 33 vents formed during Qm6 and 71 during Qm5. This indicates volcanic events occurred, on average, every 50 years during Qm5 and every 30 years during Qm6. Combined, the recurrence rate is approximately one event every 43–50 years, provided that each vent represents a separate event. Many vents are part of complex structures, often erupting simultaneously along fissures that can be up to 3 km long. For example, six closely spaced vents in the youngest Qm6 phase formed along a NW-SE fissure matching recent earthquake patterns. If an average of five vents defines one eruption episode, the true recurrence doubles to roughly every 100 years: with ten vents per event, it’s about every 200 years. These rates have remained steady since Qm3, despite uncertainties in older eruption timing.

Recurrence rates for eruptions vary from 50 to 4,000 years, depending on the period examined. These numbers include only actual eruptions, excluding failed ones, which are difficult to track due to the limited availability of historical seismic data. Since events like the 2009 failed eruption may influence future risk, it is wise to include them in probabilistic models despite data gaps. Using conservative estimates, if the 2009 event is incorporated into a model starting from 2024 and assuming a 50-year recurrence interval, an eruption could occur within the next 35 years. This should be regarded as a worst-case scenario. Locating future events in a monogenetic volcanic field, such as Harrat Lunayyir, is challenging because new eruptions are likely to occur outside known vent areas but still within the overall zone. As volcanic fields develop and more vents form, future eruptions may break through existing volcanoes, forming complex mashups — a pattern seen in western Saudi Arabia's harrats, where younger volcanoes sometimes overlap older ones with big gaps in between.

After conducting both a desktop study and an on-site reconnaissance mapping of volcanic features throughout the Lunayyir Volcanic Field, we have found that a volcanic eruption would likely unfold through a specific sequence of events observed at most of the studied locations. These eruptive processes generally follow a predictable pattern, resulting in a series of cascading volcanic hazards that should be included in any hazard assessment. We identified this chain of events and ranked its likelihood using a straightforward four-tier system. Table 1 presents the potential progression of cascading volcanic hazards along with estimated likelihoods. It's important to note that this analysis is based on rapid assessment and does not reflect a detailed or quantitative evaluation involving precise counts of cases, eruptive volumes, or durations within individual eruption periods.

In Table 2 the identified volcanic hazard estimates are demonstrated within the target area.

5. Conclusions

Our research indicates that development within the Option 4 area would significantly compromise the region’s geoheritage integrity, as it encompasses one of the most pristine environments in Harrat Lunayyir. Such activity could lead to loss of geodiversity and decrease the value of geoheritage assets. The area is especially sensitive from a geoheritage perspective [87], characterized by ash deposits and its proximity to an exceptionally preserved scoria cone. Even with careful planning, development carries inherent risks to these unique features. The expert team from the Saudi Geological Survey (SGS) has determined that the combination of moderate volcano-seismic hazard intensity—posing significant challenges for infrastructure—and the outstanding geoheritage of the site justifies a cautious approach to any further development of Option 4. Although volcanic and seismic threats are rated as low to medium in intensity, they still pose significant risks, including potential economic impacts that have already been identified. SGS’s assessment reveals that a significant portion of the Option 4 site is situated within zones unsuitable for safe development due to persistent, high-impact hazards and the absence of feasible methods to delineate secure areas.

Exclusion zones have been established in areas with high volcanic seismic activity. SGS recommends placing development projects outside these sensitive areas. Multi-night stays should be avoided; instead, controlled daytime visits may be allowed only when volcano-seismic conditions are stable. A comprehensive mitigation plan (Table 3) will help reduce potential economic losses from geological hazards and enable quick interventions during visits.