Spatio-temporal forecasting of future volcanism at Harrat Khaybar, Saudi Arabia

Abstract

The 180,000 km² of Arabian lava fields (“harrats” in Arabic) form one of the largest distributed basaltic provinces in the world. The most recent eruption in 1256 AD, on the outskirts of Medina, as well as shallow dike emplacement in 2009, ~ 200 km northeast of the city, suggest future volcanic threat to this area. Harrat Khaybar (~ 1.7 Ma to present) is one of the largest and most compositionally diverse Arabian lava fields; it is located ~ 137 km northeast of Medina and covers ~ 14,000 km². Here, we present a new eruption event record and the first estimation of future potential locations and timing of volcanism in Harrat Khaybar. Volcanic vents and eruptive fissures were mapped using remote sensing and field studies, and categorized into a geospatial database, complemented by 16 new ⁴⁰Ar/³⁹Ar ages. Our analysis reveals that Harrat Khaybar developed over five eruptive phases, where vent locations over time focus towards the central axis forming a broad N-S trend, with a central group concentrated along an axis of the regional Makkah-Madinah-Nafud (MMN) line and wider spatial dispersion between vents outwards from there. For the whole field, we estimate a long-term average recurrence rate of ~ 2.3 eruptions per 10 kyr assuming a Poisson distribution for inter-event times, which indicates that Harrat Khaybar would belong to a global group of highly active distributed volcanic fields. Our analysis also reveals that the field likely had a “flare-up” period between 450 and 300 ka where the vast majority of eruptions occurred, with ~ 18 eruptions per 10 kyr. After this intense period, eruption rates fell to < 2 eruptions per 10 kyr. Based on our findings, we estimate cumulative probabilities of 1.09 and 16.3% as lower and upper bounds of at least one eruption occurring over the next 100 years somewhere in Harrat Khaybar, with the highest probabilities within the central axis region, in particular around Jabal Qidr, Bayda and Abyad.

Introduction

The highest density of subaerial volcanic vents on Earth are found in distributed volcanic fields (basaltic volcanic fields; Valentine and Connor 2015). These volcanic fields are usually referred to as “monogenetic”, as every new eruption typically occurs at a new location. Distributed volcanic fields can be found in every tectonic environment, and many active or potentially active fields are within or near large population centers such as Auckland Volcanic Field (Auckland, New Zealand), Newer Volcanics Province (Melbourne, Australia), Boring Volcanic Field (Portland, Oregon), and Chichinautzin volcanic field (Mexico City, Mexico) (Evarts et al. 2009; Nieto-Torres and Martin Del Pozzo 2019; Heath et al. 2020; Hopkins et al. 2020). Given their longevity, variable eruptive style, and uncertain spatial patterns, understanding the volcanic hazard associated with distributed volcanic fields is an ongoing challenge (e.g., Bebbington and Cronin 2011; Connor et al. 2012; El Difrawy et al. 2013; Becerril et al. 2014, 2017; Bertin et al. 2019), which relies upon comprehensive estimations of both the location and the timing of the next eruptive event. In consequence, central to these efforts is the need to develop a quantitative understanding of the spatial and temporal aspects of the volcanic and eruptive history of a distributed volcanic field.

A common approach to determine the most likely areas of future vent opening is to convert the spatial volcano-structural data of eruptive vents (e.g., vent locations, faults, eruptive fissures, dykes) into probability density functions (e.g., Becerril et al. 2013; El Difrawy et al. 2013; Connor et al. 2019; Bertin et al. 2019), whereas the behavior and frequency of future volcanic events can be forecast by analyzing the eruption record (Bebbington and Cronin 2011; Runge et al. 2014). All this means that the spatial and temporal components of hazard are assumed as independent, so they are evaluated separately and then combined into a single model (e.g., Bebbington and Cronin 2011; Bebbington 2013; Connor et al. 2013; Bertin et al. 2019), although we acknowledge that more sophisticated approaches have tried to explicitly include the temporal component into the spatial distribution of eruptions (e.g., Connor and Hill 1995; Cappello et al. 2013). Here, to keep exploring this approach, we investigate the spatial and temporal patterns of volcanism of Harrat Khaybar, Saudi Arabia, one of the several distributed volcanic fields that make up the Arabian lava fields (Fig. 1). Despite the recent activity, the volcanic history of the harrats remains poorly known, especially for those able to impact population centers such as Rahat, Lunayyir and Khaybar. Harrat Khaybar, one of the least studied fields, is one of the most extensive volcanic fields in the Arabian Peninsula covering ~ 14,000 km². Harrat Khaybar is also characterized for producing a very wide range of eruption compositions and styles, from long lava flows (up to 58 km-long) to explosive pyroclastic density current-producing events (Camp et al. 1987).

Fig. 1

Location map of Arabian harrats in western Arabia (after Coleman et al. 1983; Camp and Roobol 1992). The gray region shows the Precambrian Arabian shield in Saudi Arabia while yellow areas represent the Cenozoic volcanic fields. The MMN volcanic fields are highlighted in orange. The green box highlights the location of the coalesced harrats: Harrat Kura, Khaybar and Ithnayn, and indicates the location of Fig. 2a. The Afar mantle plume is shown in a red circle on western Arabia (after Chang and Van der Lee 2011). Solid pink lines show plate boundaries and major tectonic features (after Stern and Johnson 2010). The blue star on the northernmost part of Harrat Rahat indicates the location of the 1256 AD eruption

In this study, we utilize remote sensing and Geographic Information System (GIS) techniques to categorize and analyze surficial volcanic characteristics (e.g., vent location, structure, morphology, alignment, density) complemented by field observations and ⁴⁰Ar/³⁹Ar geochronology. From these results, we derive the first eruption event record for Harrat Khaybar, as an initial contribution to better forecast future volcanic eruptions in the region.

Background

The Arabian lava fields, also known as harrats, form one of the largest distributed basaltic provinces in the world (Camp and Roobol 1989). These extensive volcanic fields span the western side of the Arabian Peninsula from Yemen in the south through Saudi Arabia and Jordan to Syria and Turkey in the north over ~ 180,000 km² (Fig. 1) (Camp et al. 1991). Most of these volcanic fields are late Cenozoic in age (< 12 Ma) and overlie the Precambrian Arabian shield. At least 21 eruptions took place during the last millennium on the Arabian Peninsula, based on historical records, the most recent apparently in 1937 at Dhamar in North Yemen (Siebert et al. 2011).

Saudi Arabia hosts 13 separate lava fields covering a total area of ~ 90,000 km² (Fig. 1) (Coleman et al. 1983). These are all dominated by alkali basalt, but each has an independent volcanic history and can display a variety of volcanic styles from dominantly effusive to rare local explosive (Strombolian, Hawaiian, and phreatomagmatic) activity. The most recent eruption in Saudi Arabia was in 1256 AD on the northern part of Harrat Rahat on the outskirts of Medina (Fig. 1) (Camp et al. 1987; Lindsay and Moufti 2014), whereas seismicity and ground deformation in 2009 was associated with shallow dyke emplacement in Harrat Lunayyir near Umluj on the western coast of Saudi Arabia (Pallister et al. 2010). These historical events, along with records of other late Holocene eruptions, indicate that several of the Saudi Arabian harrats are still active, and highlight the importance of conducting hazard analysis and risk assessments for these harrats to reduce the potential risks of future eruptions (Lindsay and Moufti 2014).

Regional geology

Magmatism in the Western Arabian Peninsula is temporally divided into two prominent main phases showing different chemical compositions and structural settings (Camp and Roobol 1992). The first and older period (ca. 30 to 20 Ma) was associated with passive mantle upwelling, contemporaneous with the opening of the Red Sea (Fig. 1). This produced eruptions of tholeiitic to transitional lavas along a structural trend that lies parallel to the Red Sea (Fig. 1) (Camp and Roobol 1992). A subsequent period of volcanism (ca. 12 Ma to Present) is attributed to the Afar mantle plume beneath western Arabia, developing extensional stresses and asthenospheric uplift (Fig. 1) (Chang and Van der Lee 2011; Konrad et al. 2016; Lim et al. 2020). This second and younger phase is more alkalic, producing basanite to alkali olivine–basalt and rare more evolved magmas, erupted along a N-S axis (Camp and Roobol 1992).

The largest volcanic fields of the post 12 Ma stage include harrats Rahat, Khaybar, and Ithnayn, which form a northward-younging trend (Camp and Roobol 1989; Camp et al. 1991). The central parts of the harrats are higher in elevation due to stacked lavas and scoria cones that form a prominent, N-S 600 km-long linear vent system that is termed Makkah-Madinah-Nafud (MMN) (Fig. 1) (Camp et al., 1991; Camp and Roobol 1992). The MMN vent system hosts magmas with high magmatic differentiation and a number of volcanic landforms such as the explosive volcano, Jabal Abyad, in central Harrat Khaybar. Central Harrat Khaybar erupted a suite of rocks varying from alkali-basalt to comenditic rhyolite, and Harrat Rahat erupted alkali-basalts to mugearites and trachytes (Murcia et al. 2017). The harrats flanking the MMN axis are smaller and younger and show less magmatic diversity. For instance, Harrat Lunayyir (ca. 600 ka to present) is basanite to trachy-basalt in composition, and Harrat Hutaymah (ca. 850 ka to present) is alkali basalt to tephrite (Duncan and Al-Amri 2013; Duncan et al. 2016) (Table 1).

Table 1 Comparison between several Arabian Harrats in western Saudi Arabia

Harrat Khaybar

Harrat Khaybar is a basalt-dominated volcanic field located ~ 137 km NE of Medina city, covering ~ 14,000 km² on the western Arabian Peninsula. Camp et al. (1991) assigned an age range of ca. 5 Ma to present for the exposed volcanic deposits in Harrat Khaybar based on 33 K-Ar ages.

Harrat Khaybar borders and spatially overlaps harrats Ithnayn and Kura (Camp et al. 1991) (Fig. 2a). These three distributed fields differ slightly in chemical composition. Harrat Khaybar and Ithnayn basalts are mainly mildly alkalic, producing olivine transitional basalts (OTB), whereas Harrat Kura is more alkalic, producing alkali olivine basalts (AOB) and hawaiites (Camp et al. 1991).

Fig. 2

Stratigraphic map and column of Harrat Kura, Khaybar, and Ithnayn (modified from Camp et al. 1991). a Map shows the basalt-distribution of the integrated harrats, and the dashed black box indicates the location of Fig. 3. Location of this region is illustrated in Fig. 1, the green box. b Column depicts the distribution of stratigraphic units in each basaltic field with K-Ar ages, and the red rectangle indicates Harrat Khaybar’s stratigraphic units

Collectively, harrats Khaybar, Ithnayn, and Kura, are divided into four stratigraphic units, based on disconformities, K-Ar dating and surface morphology. These include: Kura Basalt (ca. 10–5 Ma), Jarad Basalt (ca. 5–3 Ma), Mukrash Basalt (ca. 3–1 Ma), and Abyad Basalt (< ca. 1 Ma) (Table 1: Camp et al. 1991). While many of the fields contain only one or two of these units, Harrat Khaybar is a distinctive and complex distributed field that spans three stratigraphic units (Fig. 2b): Jarad, Mukrash, and Abyad basalts (Camp et al. 1991).

Harrat Khaybar is characterized by a wider variety of compositions and eruptive styles than any of the other Harrats. Camp et al. (1991) suggested that volcanic activity in Harrat Khaybar started primarily in the center of the field along a major N-S fissure producing Jarad and Mukrash OTB and forming the extensive and voluminous “whaleback” lava flows that are mainly directed to the west. Collapse craters are abundant above the extensive lava tubes of Jarad and Mukrash basalts, which resemble small volcanic vents (i.e., shield-like extrusions) (Camp et al. 1991). Then, the more differentiated rocks of Abyad basalt are attributed to a high-level magma chamber below the central spine of Khaybar that formed after the major open fissure was sealed by the first two basaltic units (i.e., Jarad and Mukrash). The central axis of the field is now defined by a N-S aligned chain of scoria cones and trachyte domes, each ~ 100 m high. Shield volcanoes, each ~ 150 m-high, are spread throughout the harrat (Camp et al. 1991). In addition, there are several felsic domes within the central area (e.g., Jabal Abyad) along with various tuff cones (e.g., Jabal Bayda) and a steep mafic composite cone (Jabal Qidr) (Fig. 3). Vents are distributed in a 100-km-long linear vent system, along a central axis with more evolved volcanics towards the center of the vent system.

Fig. 3

Aerial photograph of central Harrat Khaybar (image obtained from Google Earth). It shows the three distinctive volcanic vents resulting from different eruptive styles; Jabal Qidr (composite cone), Jabal Abyad (felsic dome), and Jabal Bayda (tuff cone) as well as other vents and lava flows in the center of the field. Location of this region is illustrated in Fig. 2a, the dashed black box

The most distinctive and best-known landforms in the region are Jabal Abyad, Bayda, and Qidr (Fig. 3) (Camp et al. 1991). Jabal Abyad dome has a diameter of ~ 600 m and reaches 2093 m a.s.l., the highest elevation amongst all the surrounding harrats. Abyad produced a thick, short coulée to the west, comprising flow-banded layers of white comendite and gray-black obsidian (Baker et al. 1973; Camp et al. 1991). Jabal Bayda, located ~ 2.5 km west of Jabal Abyad is a tuff ring comprised of poorly sorted lapilli-rich pyroclastic density current deposits, with a ~ 1500-m in diameter circular base and a small vesicular comendite dome in the center of the crater floor (Camp et al. 1991). Jabal Qidr is distinctive in this region (and Saudi Arabia) as a composite cone (~ 2022 m a.s.l.), with a ~ 400 m diameter summit crater (Camp et al. 1991). It produced ~ 18 km-long basaltic pahoehoe lava flows, which mainly flowed southwest, and a large fan of fallout lapilli that accumulated > 10 cm of scoria up to 15 km east of the summit. These three volcanic landforms in the center of the harrat represent the youngest known phases of volcanism in Harrat Khaybar (Camp et al. 1991).

Harrat Khaybar lacks detailed analysis of its volcanic history upon which estimates of its future hazard potential might be based. Along with volcanic styles and erupted volumes, it is important to identify any spatio-temporal eruption patterns, such as alignments, migration of activity, and flare-ups (e.g., Von Veh and Németh 2009; Cebriá et al. 2011; Richardson et al. 2013). Past published age data from Harrat Khaybar (Fig. 2b) were determined using K-Ar methods. Recent studies of nearby harrats Lunayyir and Hutaymah, using higher resolution ⁴⁰Ar/³⁹Ar incremental heating methods, show much younger ages than the K-Ar determinations (Duncan and Al-Amri 2013; Duncan et al. 2016). Thus, we used this new age determination technique to refine our understanding of the temporal development of Harrat Khaybar.

Methods

To determine the eruption event record for Harrat Khaybar, we carried out a morphological analysis using satellite imagery at various scales to integrate the new radiometric ages, vent locations, and tectonic structures, information that was later used to determine the relative eruption timing, volcanic alignments and spatial density estimates. The morphological analysis was conducted following the Geographic Information System methodological approaches of Kereszturi et al. (2016), Haag et al. (2019), and Morfulis et al. (2020), and integrated into a geospatial volcanic database. This database (available as Additional file 1) forms the basis of our analysis to reveal the spatial and temporal evolution of the harrat (cf. Bebbington and Cronin 2011; El Difrawy et al. 2013; Bertin et al. 2019).

Identifying volcanic vents and fissures

We created a new digital spatial database for Harrat Khaybar by identifying possible eruptive vents and fissures through high-resolution satellite images (~ 30 m) based on geometric and visual criteria using ArcGIS Pro software (e.g., Haag et al. 2019). Besides satellite images, we used a digital elevation model (DEM) generated from the 30-m Shuttle Radar Topography Mission (SRTM; available at https://earthexplorer.usgs.gov/) to derive topographic contours, slope, aspect, and hillshade maps (cf., Haag et al. 2019). We compared our results with the 1:250,000-scale geologic maps of Roobol and Camp (1991) and a previous catalogue of vents from Runge et al. (2016). Finally, we analyzed rasterized, reduced-to-pole aeromagnetic data of Zahran et al. (2003) to investigate the presence of buried vents at Harrat Khaybar (cf., Yucca Mountain area, O'Leary et al. 2002; Harrat Rahat, El Difrawy et al. 2013).

After mapping volcanic vents and fissures in the harrat, we catalogued them into temporal eruptive phases. We incorporated the new radiometric data into our new digital spatial database and then looked for relations between the mapped volcanic vents in terms of morphological preservation and superpositional relationships. Vents that showed similarities in these criteria were then assigned to one specific group.

Age data

We collected samples in two field campaigns for age determinations (Fig. 4). Samples 1–8 are from an east-west transect aimed at dating the full stratigraphy of Jarad, Mukrash, and Abyad volcanic units, while samples 9–16 are concentrated in young, elevated volcanic vents aimed at dating the evolved (felsic) compositions. Rocks were broken from flow interiors, far removed from weathered surfaces.

Fig. 4

Aerial photograph shows the locations of samples collected for age determinations

We performed age determinations of groundmass separates at Oregon State University with standard ⁴⁰Ar/³⁹Ar laser step heating methods. Whole rock samples were crushed, sieved, washed, and subjected to mild acid leaching with HCl and HNO₃ then hand-picked to remove visible phenocryst-containing fragments. Some 100 mg of prepared groundmass of each sample was next irradiated for 6 h in the 1 MW TRIGA nuclear reactor at Oregon State University, along with flux monitor FCT sanidine (28.201 Ma; Kuiper et al., 2008). We loaded irradiated samples into Cu-planchettes in an ultra-high vacuum sample chamber and incrementally heated (in 30–40 steps, from 400 °C to fusion) by scanning a defocused 25 W CO₂ laser beam in preset patterns across each sample, in order to release the Ar evenly. Argon isotope compositions of irradiated samples were determined using the Thermo Scientific Model ARGUS VI multi-collector with five fixed Faraday detectors (all fitted with 10¹² Ω resistors) and 1 ion-counting CuBe electron multiplier, allowing simultaneous measurement of all Ar isotopes, with mass 36 on the ion multiplier and masses 37 through 40 on the four adjacent Faraday cups. Gas cleanup occurred in a small-volume, all-metal extraction line equipped with Zr–Al getters. We monitored the atmospheric correction with an air pipette system, and cross-calibrated the 5 collectors for small differences in sensitivity, on a daily basis. We calculated ages (reported at ±2σ uncertainty) in several ways. Plateau ages are the weighted (by inverse variance) means of concordant, sequential step ages while isochron ages are derived from regressions of the same heating step isotopic compositions. All calculations were performed using the ArArCALC v2.6.2 software package (Koppers, 2002). Full datasets, including spectra (plateau) and inverse isochron plots, are reported in Additional file 2.

In all sample analyses (but HK-5), plateau and isochron ages are concordant, indicating simple Ar release from essentially fresh groundmass during step heating. In many samples, the initial (trapped) ⁴⁰Ar/³⁶Ar calculated from isochrons is of atmospheric composition (298). In some samples, the calculated initial ⁴⁰Ar/³⁶Ar is slightly less than the atmospheric value, indicating fractionation of ³⁶Ar from ⁴⁰Ar, probably into vesicles as the lavas cooled, and in a few samples, the initial ⁴⁰Ar/³⁶Ar is slightly greater than the atmospheric value, indicating ‘excess’ or undegassed, mantle-derived Ar. We calculate precise plateau ages from consecutive age-equivalent steps, comprising 50–100% of total gas released. MSWD, the mean square of weighted deviation, is an F-statistic that compares variation within step ages with variation between step ages, and all measures (except for HK-5) indicate significance of the plateau ages. Sample HK-5, a comendite, produced a disturbed age release pattern (no statistically significant plateau) and is not reported in the Additional file 2.

Alignment analysis

Detecting alignments in volcanic fields using point features (i.e., volcanic vent coordinates) can shed some light on the key controls of intraplate volcanic origin or migration of magma. In other words, vent alignments can represent areas of magma focusing or structural weaknesses, along which volcanism returns in repeated episodes (e.g., Kear 1964; Connor et al., 2000; Le Corvec et al. 2013; Morfulis et al. 2020; Sieron et al. 2021). Several methods have been developed to identify orientations in volcanic fields, including: two-point azimuth techniques (Lutz 1986; Lutz and Gutmann 1995; Bleacher et al. 2009; Richardson et al. 2013), spatial transform methods (e.g., Hough’s transform; Connor 1990; Connor et al. 1992; Von Veh and Németh 2009), and strip methods (Zhang and Lutz 1989; Arcasoy et al. 2004). In our study, we focused on the nearest neighbor azimuth method of Cebriá et al. (2011) using the Thomson and Lang (2016) Matlab tool (e.g., Grosse et al. 2020; Morfulis et al. 2020). This technique is an improved version of the two-point azimuth method of Lutz (1986). It is based on the concept of nearest neighbor following the ideas of Kear (1964), which suggest that any monogenetic volcanic vent is more likely to be generated by the same fissure that fed its nearest vent. Thus, straight lines between two neighboring vents within a determined distance (i.e., line length) would represent meaningful lineaments. The determined distance is optimally obtained using the concept of the minimum significant distance. For local scale features, it is suggested to be equal to the third lower part of the anomalous end of all distances between all vents (i.e., d ≤ ([x̄- σ]) /3), where d is the selected distance, σ is standard deviation, and x̄ is the mean of distances (Cebriá et al. 2011). For regional-scale structures, more widely spaced vents should be considered. The frequencies of azimuths of vent-connecting lines can be displayed in a rose diagram, where azimuths that have a frequency higher than one standard deviation above the mean (i.e., x̄ + σ) are more likely to represent the preferred orientations of vent alignments. Furthermore, the modified two-point azimuth method of Cebriá et al. (2011) can also recognize curved features defined by the combination of several short straight lines.

Spatial probability analysis

After mapping volcanic vents and describing eruptive stages, we performed a spatial probability analysis to estimate the density of visible vents and forecast the distribution of future volcanic activity using a kernel density estimation (cf., Connor et al. 2015; Connor et al. 2019). This is a non-parametric statistical technique used for calculating probability density function of point-like features (Connor et al. 2019). In this study, we utilized the MatHaz program of Bertin et al. (2019). MatHaz is an open-source Matlab code to conduct probabilistic spatio-temporal volcanic hazard assessments in distributed volcanic fields. Probability isocontours in this process are calculated using Gaussian kernel smoothing (Bertin et al. 2019). Elliptical bandwidth estimators can also be selected, such as the Asymptotic Mean Squared Error (AMSE) and the Sum of Asymptotic Mean Squared Error (SAMSE) (Duong 2005). Data included here were only vent locations due to limitation of data in Harrat Khaybar. To take into account the temporal history of the field, we applied a decreasing exponential age-weighting procedure, so that youngest vent-age information is most important to the final kernel than old vents (see Bertin et al. 2019). The bandwidth matrix used in this study for all vents is:

$${\hat{\textrm{H}}}_{\textrm{PI},\textrm{AMSE}}=\left[\begin{array}{cc}5497455& -3560658\\ {}-3560658& 61844432\end{array}\right]$$

chosen using the AMSE optimization method (Duong 2005). This corresponds to an ellipse with major and minor axis lengths of 7.88 and 2.30 km, respectively, and whose major axis is rotated 4° to the west from the north (i.e., an azimuth of 356°). Eventually, the spatial probability analysis is reported as a matrix, where each cell has a probability.

Temporal probability analysis

To determine temporal recurrence rates for the eruptive stages and overall volcanism at Harrat Khaybar, we applied both a long-term average recurrence rate and power law process using MatHaz. The long-term average recurrence rate follows a simple Poisson distribution for the repose-period between eruptions (e.g., Connor et al. 2013; El Difrawy et al. 2013; Runge et al. 2014). This is based on the assumptions that the recurrence rate is constant, and that each vent represents a separate eruption (e.g., Ho et al. 1991). It can be obtained by knowing the number of volcanic vents/eruptions over a specified time interval (Ho et al. 1991; Bertin et al. 2019), which is written:

$$\phi\ \left(\Delta t\right)=\frac{N-1}{t_1-{t}_N}$$

Where N is the total number of eruptive vents that falls between t₁ (the age of the oldest known vent) and t_N (the age of the youngest known vent).

The second approach used is the power law process or Weibull distribution. It follows a nonhomogeneous Poisson distribution that accounts for changes in eruption rates, commonly used to describe accelerated failure with time (Connor et al. 2015). This approach requires definite ages for all eruptive vents (Ho et al. 1991). One of its drawbacks is that it might produce physically unreasonable results, especially at longer time intervals (Connor et al. 2015).

For the long-term average recurrence rate approach, we assigned a relative age range to each age group (i.e., eruptive episode) based on the available absolute ages. For the power intensity function approach, we estimated definite ages (with no age uncertainties) for each vent based on the identified age ranges. This was possible by generating random variates from a discrete uniform distribution with a seed number of 999 using Crystal Ball in Excel (e.g., Connor et al. 2013).

Spatio-temporal analysis

Due to the lack of a robust geochronological database for Harrat Khaybar, we must assume that both the spatial and temporal components of volcanism are independent (cf., Bebbington and Cronin 2011; Bebbington 2013; Connor et al. 2013; Bertin et al. 2019). Each component is calculated separately and then multiplied by each other to obtain the spatio-temporal estimation of future vent opening. In other words, we multiplied the matrix (spatial analysis) by a scalar (the temporal analysis/recurrence rate). The probabilities of occurrence are then obtained by integrating the spatio-temporal matrix over a finite region, which in our case was our whole study area.

Results

Geochronology

⁴⁰Ar/³⁹Ar data obtained for 16 samples are summarized in Table 2. Complete data files for all experiments are available in Additional file 2. Ages determined to bracket the mapped stratigraphic units are from seven alkali basalt lava flows and one comendite (HK- samples), while the other eight ages are from rocks with more felsic compositions within the Abyad unit (K- samples). Reliable plateau ages for 15 of these units range from 1.7 Ma to 53 ka. The majority of the 31–40 incremental heating steps analyzed in each of the experiments (50–100% of total ⁴⁹Ar released) define the age plateaus, and the MSWDs comparing variance within step ages versus variance of step ages about the plateau age are acceptably low for significance of plateau ages. These have corresponding concordant isochron ages with atmospheric ⁴⁰Ar/³⁶Ar intercepts. One sample (HK-5 Mukrash comendite), however, displayed an erratic step heating release pattern resulting in a very large MSWD (> 100); based on this we do not consider its age as reliable as the others. We see no evidence for excess ⁴⁰Ar, as reported in the age spectra for several samples from nearby Harrat Hutaymah, which has been linked to mantle and crustal xenolithic fragments and the likely cause for K-Ar ages being generally older than more recent ⁴⁰Ar/³⁹Ar incremental heating ages (Duncan et al. 2016).

Table 2 Age determinations for volcanic rocks from Harrat Khaybar, western Saudi Arabia. Map# refers to circled numbers on Fig. 4

The new data from the stratigraphic transect indicate the Jarad activity began at or before 1.7 Ma, followed by Mukrash activity between 1.5 and 0.35 Ma, and then by Abyad activity from 0.35 Ma to present, where central volcanoes developed felsic compositions from at least 0.95 Ma to present. Although Camp et al. (1991) define their stratigraphic units partly on disconformable erosional surfaces, we see no obvious age breaks in volcanic activity. Our ⁴⁰Ar/³⁹Ar age ranges are distinctly younger and tighter than the K-Ar ages reported by Camp et al. (1991) (Jarad, 5–3 Ma; Mukrash, 3–1 Ma; Abyad 1 Ma to present). This has important implications for revising the frequency of volcanic eruptions.

Eruption event record

Three hundred eighty-four volcanic vents were identified in our study at Harrat Khaybar. Whereas volcanic events might be associated with multiple vents (e.g., Condit and Connor 1996; Muffler et al. 2011; Runge et al. 2014; Gallant et al. 2018, 2021; Nieto-Torres and Martin Del Pozzo 2019; Downs et al. 2020), in this work, a vent is assumed to represent an individual volcanic event (e.g., Bebbington and Cronin 2011; El Difrawy et al. 2013; Bertin et al. 2019; Connor et al. 2019). These vents belong to different volcanic morphologies, including: scoria cones, small-shield shield volcanoes, tuff cones, and a composite volcano. Of these vents, 83 define 37 possible eruptive fissures. All spatial data are available in Supplementary files S1A and S1B in Additional file 1.

Our new age data indicates an age range of ca. 1.7 Ma to present for Harrat Khaybar (Table 2), where about 84% of the mapped volcanic vents erupted < 600 ka. We divide this history into five age groups (Age-1 to Age-5 here, with Age-1 being the oldest group) based on identifying and grouping vents by morphological preservation and superpositional relationships, and then establishing best-fit age ranges for each group (Table 3, Fig. 5). The pre 600 ka vents are highly eroded and show a wide age range while the post 600 ka vents are relatively well-preserved and have better temporal constraints. Therefore, we break the last 600 ka eruptions into four even size periods although we acknowledge that there are significant uncertainties in these eruptive periods since only ~ 4% of the vents have definitive ages, which is a common situation in poorly dated volcanic fields (e.g., Connor et al. 2013; El Difrawy et al. 2013; Gallant et al. 2018, 2021; Nieto-Torres and Martin Del Pozzo 2019; Valentine et al. 2021).

Fig. 5

Spatial distribution of volcanic vents at Harrat Khaybar. a, b, c, d, and e Maps represent the spatial distribution of Age-1, 2, 3, 4, and 5 groups, respectively. f Map shows the spatial distribution of the eruptive fissures at the harrat

Table 3 Physical characteristics of volcanoes in each age group at Harrat Khaybar

Age-1 group defines the first eruptive stage that occurred between 1.7 Ma and 600 ka (Fig. 5). It consists of 62 basaltic small-shield volcanoes (i.e., 62 vent/eruptions) and three eruptive fissures. Age-2 group (600 to 450 ka) comprises 20 trachyte, benmoreite, and mugearite lava domes. Age-3 group (450 to 300 ka) includes 268 basaltic scoria cones and 28 fissures. Age-4 group (300 to 150 ka) shows only five felsic explosive eruptions, tuff cones and lava domes of comendite in the center of the harrat. The latest eruptive stage is Age-5 group (< 150 ka), which consists of 29 small-shield volcanoes, six fissures and the composite cone of Qidr (Table 3, Fig. 5).

Along with the temporal frequency variation, the spatial distribution of volcanic vents of each eruptive phase varies over time as well. Volcanism apparently began with scattered shield volcanoes across a broad area, then it concentrated in the central field from the second stage and focused on increasingly narrow areas. The third eruptive period generated a clear linear vent system, with the vast majority of vents in the whole field history concentrated along the central axis. The fourth period erupted only in the center of the harrat, while the most recent volcanic activity involved shield volcanoes and a composite cone along a short segment of the central alignment, (Figs. 4a to f, Table 3). Kernel density estimations conducted for each eruptive phase (except Age-4 Group which has too few vents for analysis) and then analyzed together, show that volcanism at Harrat Khaybar has a tendency to cluster in the central part of the field (see Additional file 3). Likewise, the spatial probability analysis conducted for all the vents (384), using an age-weighting procedure, denotes an elongated area of high probability density (between 10^− 5 and 10^− 4 vents/0.25 km²) along the center of the harrat (Fig. 6).

Fig. 6

Spatial probability isocontour map of all visible volcanic vents. It depicts the distribution of probabilities throughout the harrat based on the locations of all visible vents using an age-weighting procedure in order to give less weight to older vents

The available aeromagnetic data for Harrat Khaybar did not provide insight into the presence of buried volcanic vents, as there was no evident correlation between magnetic anomalies and locations of visible vents (see Additional file 4) nor any differentiation between the volcanics and the Precambrian basement. There was also no other supportive data for the presence of buried vents at Harrat Khaybar such as lava flows thickness and seismic profiling. At higher resolution, it might be useful, but at this scale, it only supports the N-S alignment of the volcanic field as elongated, positive magnetic anomalies along the central spine of the harrat.

Alignment analysis

The quantitative alignment analysis revealed an overall NNW orientation during the majority of the eruptive periods at Harrat Khaybar. Age-1 group shows a polymodal tendency alignment with the highest frequency towards the NW (Fig. 7a). Age-2 group is organized along predominant NW-directed vent alignments (Fig. 7b). Age-3 and 5 groups show NNW strikes (Fig. 7c-d). The fourth period has insufficient vent data for analysis. The predominant orientation of all vents is N-S (Fig. 7e, Table 3). These vent alignments obtained by the nearest neighbor azimuth method are consistent with the azimuth of major axes derived from the spatial probability analysis, except alignments of Age-1 and 2 groups that show significant differences in directions (see Supplementary file S1C, Additional file 1).

Fig. 7

Rose diagrams of azimuths obtained from vents-connecting lines. Each diagram represents the strike frequency of the mapped lineaments at 10° bin intervals. Blue bins are those with frequencies higher than one standard deviation above the mean (σ + x̄). The black arrow represents the vector mean

Temporal probability analysis

The previously determined five age ranges that span from ca. 1.7 Ma to present are used in estimating the average long-term recurrence rate of each age group (Table 3) (see Supplementary file S1D, Additional file 1). According to Fig. 8a, there is a clear fluctuation by two orders of magnitude in recurrence rates among age groups, where the peak of activity occurred during Age-3 group (450 to 300 ka) with a recurrence rate of 17.8 eruptions/10 kyr (Table 3). The long-term average recurrence rate for the whole field is 2.25 eruptions/10 kyr, which places it within a group of high-rate fields around the world using their number of vents and lifespan (e.g., Springerville in the USA (2.78 eruptions/10 kyr), Northern Harrat Rahat in Saudi Arabia (2.91 eruptions/10 kyr), Armenia (3.46 eruptions/10 kyr), and Eastern Snake River Plain in the USA (4.22 eruptions/10 kyr)) (Fig. 8b; Connor and Hill 1995; Connor and Conway 2000; Weller 2004; El Difrawy et al. 2013; Valentine and Connor 2015; Gallant et al. 2018; Valentine et al. 2021; and references therein).

Fig. 8

Charts showing eruption recurrence rates of age groups in Harrat Khaybar and average rates of selected basaltic volcanic fields around the world. a Chart of the recurrence rates among age groups at Harrat Khaybar. b Chart comparing average temporal recurrence rates of well-studied distributed volcanic fields in the world including Harrat Khaybar. All data points are averages calculated from the known number of vents divided by the longevity of the field. Chart a) demonstrates that average rates shown in b) do not adequately characterize the temporal variation of activity in a field. Data sources are from Connor and Hill (1995), Connor and Conway (2000), Weller (2004), El Difrawy et al. (2013), Valentine and Connor (2015), Gallant et al. (2018), and Valentine et al. (2021)

Spatio-temporal probability analysis

The spatio-temporal estimations of future volcanic activity for Harrat Khaybar are obtained by multiplying the values of the spatial analysis of the last 300 ka of activity (i.e., Age-4 and Age-5 groups) by the long-term average recurrence rate of these vent age groups (1.1 eruptions/10 kyr), and the recurrence rate of Age-3 group that is higher than other periods by about two orders of magnitude (17.8 eruptions/10 kyr), for normal and worst-case scenarios, respectively. The results of the Weibull process have been omitted from the spatio-temporal analysis because they would have added more uncertainty to the spatio-temporal probabilities as they are based on arbitrarily generated absolute ages and without age uncertainties (see Supplementary file S1D and S1E, Additional file 1). Subsequently, we utilized our findings to forecast the occurrence of at least one eruption anywhere in the field for the next 1, 10, 100 and 1000 years using the MatHaz tool (Bertin et al. 2019). This resulted in spatio-temporal cumulative probabilities of 1.09 and 16.3% (lower and upper limits) for an eruption somewhere in the harrat over the next century (Table 4; Fig. 9a and b, respectively). The spatio-temporal estimations for the next 1,10, and 1000 years using only data for the last 300 ka eruptions are available in Additional file 5. The spatio-temporal cumulative probabilities using a power law approach are available in Supplementary file S1F, Additional file 1.

Table 4 Spatio-temporal cumulative probabilities, considering the long-term average recurrence rate approach. They are obtained by multiplying the spatial density by the point process intensity function

Fig. 9

Spatio-temporal probability isocontour maps. a and b Maps show the probabilities of occurrence of at least one eruption in Harrat Khaybar for the next 100 years, obtained by multiplying the values of the spatial analysis of the last 300 ka eruptions by the long-term average recurrence rates of 1.1 eruptions/10 kyr and 17.8 eruptions/10 kyr, respectively

Discussion

This is the first study to quantitatively explore the spatio-temporal distribution of volcanism of Harrat Khaybar. The highest spatial vent densities are concentrated along a central alignment with a predominant N-S orientation, and most of the vents (268 out of 384) occurred in a single “flare-up” period between 450 and 300 ka. There has been a clear focusing of vent locations into the center of the field over time since ~ 600 kyr (Fig. 5b-e and Additional file 3), after a period of ~ 1100 kyr of dispersed volcanism. This focusing could relate to changes in the origin, migration of magma, development of pre-heated pathways, or intrusion into the crust (e.g., Kósik et al. 2020). Camp et al. (1991) proposed that a major N-S fissure was the main focus of activity at Harrat Khaybar, although initial activity does not follow this pattern. Volcanism at Harrat Khaybar is distributed, unlike other volcanic fields that show transitional activity between central-vent and distributed volcanism such as the Los Tuxtlas Volcanic Field (Veracruz, Mexico), the Tolbachik volcanic field (Kamchatka, Russia), and the Kirishima volcano group (Kyushu, Japan) (Miyabuchi et al. 2013; Kugaenko and Volynets 2019; Sieron et al. 2021). The more differentiated rocks of Abyad basalt are attributed to a shallow magma chamber below the central spine of Khaybar. The focusing of activity in age groups 2 and 3, along with the prominence of evolved compositions in the most recent stages, is consistent with the idea of preconditioning of the crust over time by successive intrusions, possibly leading to a persistent mid-crustal reservoir beneath central Harrat Khaybar. These characteristics appear to follow increased activity in age group 3 that might indicate an increased magmatic flux (i.e., a magmatic flare-up) and the consequent merging of dikes to form a unitary reservoir (e.g., Hardee 1982; Wu et al. 2021, 2022). With a growing volcanic load at the surface, dike capture and magmatic lensing could also generate a field-centered reservoir (Valentine 1993; Karlstrom et al. 2009). Crustal pre-conditioning and focusing of melts could, in-turn, create the conditions for differentiation and evolution of the system from early mildly alkalic, olivine transitional basalts (OTB) to evolved basalts, trachytes, and comendites. We suggest this transition to have occurred at ca. 300 ka.

Long term recurrent rates (number of vents divided by the longevity of activity) provide a useful synoptic comparison of distributed volcanic fields (Fig. 8). However, this should be used with caution since such fields experience waxing and waning cycles of activity, although their recurrence rates are often reported as long term (i.e., based on the whole history of the field). For instance, well-studied fields like the Springerville in Arizona, USA and northern Harrat Rahat in western Saudi Arabia have been shown to have considerable fluctuation in eruption rates during their lifespans (e.g., Mnich and Condit 2018; Stelten et al. 2020; Valentine et al. 2021, and references therein). At the Springerville volcanic field, the volume output was much higher in the early stages (< 2 to 1.50 Ma) compared to the latest stage (< 0.5 to 0.25 Ma) (Condit and Connor 1996; Mnich and Condit 2018). Closer to Harrat Khaybar, Stelten et al. (2020) show that northern Harrat Rahat has experienced 12 eruptive stages during the last 1.2 Ma years with significant magmatic flare-ups between 460 ka and 360 ka, and between 260 ka and 180 ka. The biggest pulses appear to be during eruptive stages 5, 6, and 9 with the highest recurrence rate of 6.5 eruptions/10 kyr during eruptive period 6 (260 to 180 ka), whereas the lowest recurrence rate is calculated as 0.14 eruptions/10 kyr during eruptive period 12 (780 to 1200 ka), over an order of magnitude less. Harrat Khaybar shows a similar variability of rates during the history of the field (Fig. 8a, Table 3). Of particular note is that the recurrence rate of the “flare-up” during Age-3 group (450 to 300 ka) is 17.8 eruptions/10 kyr, almost about two orders of magnitude higher than the lowest rate during Age-4 group (300 to 150 ka). The Harrat Khaybar flare-up therefore overlaps with the earlier of these northern Harrat Rahat upsurges during eruptive stage 9 and 8. This may suggest a regional increase in the flux of mantle-derived magmas into the crust during this period and a possible deep magmatic connection between two harrats. However, the major flare-up at northern Harrat Rahat during stages 5 and 6 appear to coincide with a lull during Age-4 group at Harrat Khaybar, indicating that a regional correlation may not exist throughout the history of these harrats.

One caution around the Harrat Khaybar “flare-up” is the assumption that each vent records an individual eruption, and that eruptions are of the same volume throughout the history of the field. It is common in many distributed volcanic fields that multiple vents are simultaneously produced during a single event (e.g., Muffler et al. 2011; Bevilacqua et al. 2017; Gallant et al. 2018, 2021; Connor et al. 2019; Downs et al. 2020). An example from the Arabian harrats is the 1256 AD event in northern Harrat Rahat that produced six scoria cones along a ~ 2.25-km-long fissure (Camp et al. 1987). Further, Runge et al. (2014) and Stelten et al. (2020) suggested the vents to events proportion for northern Harrat Rahat is > 1. That could also be true for Harrat Khaybar, but the problem is that applying a normalization factor to just period 3 seems subjective because the same over representation of vents relative to actual eruptions could be universal to all the time periods at Harrat Khaybar. We therefore maintain that period three represents a significant increase in activity consistent with the idea of a flare-up. Future probabilistic studies should assess the spatio-temporal relationships between vents in the harrat and consider modeling eruptive events to account for multiple vent eruptions (e.g., Condit and Connor 1996; Cappello et al. 2013; Runge et al. 2014; Gallant et al. 2018, 2021; Nieto-Torres and Martin Del Pozzo 2019). That would result in a better spatio-temporal estimation of future volcanic events.

Alignment analysis

The vent alignments using nearest neighbor azimuth analysis may indicate geologically meaningful lineaments of both local- and regional-scale features such as the underlying structural settings of the harrat (i.e., dominant stress regime) (e.g., Kear 1964; Connor et al., 2000; Le Corvec et al. 2013; Morfulis et al. 2020; Sieron et al. 2021). Cebriá et al. (2011) suggest that lineaments between neighboring vents within the minimum significant distance can reveal the local-scale faults and fissures. In fact, by analyzing vent alignments of each eruptive stage separately, slight shifts in orientations between the eruptive stages are evident (Fig. 6), which are consistent with the azimuths of the ellipse’s major axes for age groups 3 to 5 (ellipses obtained when applying the kernel density estimation method). This would imply, at a local scale, pre-existing structural conditions reflecting local stress-strain conditions (i.e., areas of crustal weakness) at the time of magma emplacement (Gonnermann and Taisne, 2015). For regional scale features, a much larger distance between neighboring vents should be considered (Cebriá et al. 2011). We have combined all vent data from this study and implemented the minimum significant distance approach to image regional stress orientation at Harrat Khaybar. This shows a prevalent N-S direction (Fig. 7e), subparallel to the N-S trending of the MMN line (Fig. 1) (Camp and Roobol 1992).

Spatial probability analysis

The area of high probability density (≥1.0 × 10^− 4 vents/0.25 km²), determined along the central alignment of Harrat Khaybar, depicts the densest zone of volcanic vents and areas more prone to experience volcanic activity in the future. This high probability zone includes a major road and a few settlements. Furthermore, the harrat is a source of water to the surrounding towns and agriculture, so volcanic activity could pose a wider issue (e.g., Valente et al. 2022). The large tephra apron of the Qidr edifice also suggests a potential aviation hazard (e.g., Sulpizio et al. 2012). The closest residential areas are Khaybar and Alhait towns, at ~ 65 km from the center of the harrat, but they have exceedingly low probabilities for new vent locations (~ 10^− 20) (Fig. 6). However, lava flows, pyroclastic density currents, and tephra fallout could impact these sites, particularly from any new vents on the edge of the low probability areas.

Temporal probability analysis

Due to the paucity of age data at Harrat Khaybar, we must account for uncertainty in determining an average recurrence rate, as this value depends on the number of eruptions over a specific time interval. To address this, we choose northern Harrat Rahat (also known as Harrat Al-Madinah) for this analysis as an analog since they share the same origin and show similar petrological characteristics (Camp et al. 1991). If we consider that Harrat Khaybar and northern Harrat Rahat have the same age range, ca. 1.7–0 Ma with 384 and 496 vents, respectively (e.g., El Difrawy et al. 2013), we obtain similar average recurrence rates 2.25 eruptions/10 kyr and 2.9 eruptions/10 kyr, respectively. However, when we consider only the ≤300 ka volcanic activity in both harrats, we obtain average recurrence rates of 1.1 eruptions/10 kyr for Harrat Khaybar and 3.0 eruptions/10 kyr for northern Harrat Rahat (El Difrawy et al. 2013). Another analogous volcanic field in terms of compositional bimodality and temporal recurrence rate is the Eastern Snake River Plain (ESRP), Idaho, USA (Kuntz et al. 1992; Gallant et al. 2018). Considering all visible vents erupted between 1.2 Ma and 2.1 ka in the ESRP (506 vents) results in an average recurrence rate of 4.22 eruptions/10 kyr, and considering only the ≤500 ka volcanic activity (286 vents) shows similar average recurrence rate 5.72 eruptions/10 kyr (Gallant et al. 2018). These average recurrence rates are within the same order of magnitude as the calculated recurrence rates for Harrat Khaybar and northern Harrat Rahat. The lower rate for Khaybar obtained here by changing the time to evaluate (1.1 eruptions/10 kyr) is a more reasonable estimate to use for forecasting future volcanic activity, as it can be more indicative of the current hazard.

Conclusion

This research aimed to determine an eruption event record and a spatio-temporal evaluation of future volcanic activity for Harrat Khaybar using remote sensing, GIS techniques, field data, and ⁴⁰Ar/³⁹Ar geochronology. Based on a quantitative and qualitative analysis of preserved volcanic vents, we conclude that Harrat Khaybar has formed during five phases of volcanism starting ca. 1.7 Ma years ago, producing at least 384 volcanic vents. The overall spatial distribution indicates a north-south trending vent system, assembled in a central alignment, that becomes more scattered towards the flanks. The spatial analysis also indicates that those areas most likely to experience volcanic vent opening in the future are found along the center of the harrat, following the identified N-S pattern.

The average recurrence rate of 1.1 eruptions/10 kyr calculated for the last 300 ka history of Harrat Khaybar provides a preliminary hazard estimate. When this information on recurrence is combined with the spatial distribution of volcanism, we obtained probabilities between 1.09 and 16.3% of at least one eruption anywhere in the field, with the central part of the harrat more likely to experience new volcanic activity compared to the flanking regions. An extreme hazard, such as a renewed flare-up in the field could see recurrence rates as high as 17.8 eruptions/10 kyr. We acknowledge that the uncertainties in current analysis are high and emphasize the need for more detailed geochronology to inform the spatio-temporal analysis.