Skip to main content

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


The 180,000 km2 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 km2. 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 40Ar/39Ar 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.


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 km2. 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
figure 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 40Ar/39Ar 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.


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 km2 (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 km2 (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 km2 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
figure 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
figure 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 40Ar/39Ar 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.


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 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
figure 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 40Ar/39Ar laser step heating methods. Whole rock samples were crushed, sieved, washed, and subjected to mild acid leaching with HCl and HNO3 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 CO2 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 1012 Ω 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) 40Ar/36Ar calculated from isochrons is of atmospheric composition (298). In some samples, the calculated initial 40Ar/36Ar is slightly less than the atmospheric value, indicating fractionation of 36Ar from 40Ar, probably into vesicles as the lavas cooled, and in a few samples, the initial 40Ar/36Ar 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 t1 (the age of the oldest known vent) and tN (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.



40Ar/39Ar 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 49Ar 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 40Ar/36Ar 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 40Ar, 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 40Ar/39Ar 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 40Ar/39Ar 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
figure 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 km2) along the center of the harrat (Fig. 6).

Fig. 6
figure 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
figure 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
figure 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
figure 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


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 km2), 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.


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 40Ar/39Ar 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.

Availability of data and materials

The datasets supporting the conclusions of this article are available in Additional files 1 (multiple files), 2, 3, 4 and 5.



Geographic Information System


Digital Elevation Model


Shuttle Radar Topography Mission


Saudi Geological Survey


  • Arcasoy A, Toprak V, Kaymakçı N (2004) Comprehensive strip based lineament detection method (COSBALID) from point-like features: a GIS approach. Comput Geosci 30(1):45–57

    Article  Google Scholar 

  • Baker PE, Brosset R, Gass IG, Neary CR (1973) Jebel al Abyad: a recent alkalic volcanic complex in western Saudi Arabia. Lithos 6(3):291–313

    Article  Google Scholar 

  • Bebbington MS (2013) Assessing spatio-temporal eruption forecasts in a monogenetic volcanic field. J Volcanol Geotherm Res 252:14–28

    Article  Google Scholar 

  • Bebbington MS, Cronin SJ (2011) Spatio-temporal hazard estimation in the Auckland volcanic field, New Zealand, with a new event-order model. Bull Volcanol 73:55–72

    Article  Google Scholar 

  • Becerril L, Bartolini S, Sobradelo R, Martí J, Morales JM, Galindo I (2014) Long-term volcanic hazard assessment on El Hierro (Canary Islands). Nat Hazards Earth Syst Sci 14(7):1853–1870

    Article  Google Scholar 

  • Becerril L, Cappello A, Galindo I, Neri M, Del Negro C (2013) Spatial probability distribution of future volcanic eruptions at El Hierro Island (Canary Islands, Spain). J Volcanol Geotherm Res 257:21–30

    Article  Google Scholar 

  • Becerril L, Martí J, Bartolini S, Geyer A (2017) Assessing qualitative long-term volcanic hazards at Lanzarote Island (Canary Islands). Nat Hazards Earth Syst Sci 17(7):1145–1157

    Article  Google Scholar 

  • Bertin D, Lindsay JM, Becerril L, Cronin SJ, Bertin LJ (2019) MatHaz: a Matlab code to assist with probabilistic spatio-temporal volcanic hazard assessment in distributed volcanic fields. J Appl Volcanol 8:4

    Article  Google Scholar 

  • Bevilacqua A, Bursik M, Patra A, Pitman EB, Till R (2017) Bayesian construction of a long-term vent opening probability map in the Long Valley volcanic region (CA, USA). Stat Volcanol 3(1):1

    Article  Google Scholar 

  • Bleacher JE, Glaze LS, Greeley R, Hauber E, Baloga SM, Sakimoto SEH, Williams DA, Glotch TD (2009) Spatial and alignment analyses for a field of small volcanic vents south of Pavonis Mons and implications for the Tharsis province, Mars. J Volcanol Geotherm Res 185(1-2):96–102

    Article  Google Scholar 

  • Camp VE, Hooper PR, Roobol MJ, White DL (1987) The Madinah eruption, Saudi Arabia: magma mixing and simultaneous extrusion of three basaltic chemical types. Bull Volcanol 49(2):489–508

    Article  Google Scholar 

  • Camp VE, Roobol MJ (1989) The Arabian continental alkali basalt province: part I. evolution of Harrat Rahat, Kingdom of Saudi Arabia. Geol Soc Am Bull 101:71–95

    Article  Google Scholar 

  • Camp VE, Roobol MJ (1992) Upwelling asthenosphere beneath western Arabia and its regional implications. J Geophys Res 97(B11):15255–15271

    Article  Google Scholar 

  • Camp VE, Roobol MJ, Hooper PR (1991) The Arabian continental alkali basalt province: part II. Evolution of Harrats Khaybar, Ithnayn, and Kura, Kingdom of Saudi Arabia. Geol Soc Am Bull 103(3):363–391

    Article  Google Scholar 

  • Cappello A, Bilotta G, Neri M, Del Negro C (2013) Probabilistic modeling of future volcanic eruptions at Mount Etna. J Geophys Res 118(5):1925–1935

    Article  Google Scholar 

  • Cebriá JM, Martín-Escorza C, López-Ruiz J, Morán-Zenteno DJ, Martiny BM (2011) Numerical recognition of alignments in monogenetic volcanic areas: examples from the Michoacán-Guanajuato volcanic field in Mexico and Calatrava in Spain. J Volcanol Geotherm Res 201(1–4):73–82

    Article  Google Scholar 

  • Chang SJ, Van der Lee S (2011) Mantle plumes and associated flow beneath Arabia and East Africa. Earth Planet Sci Lett 302(3–4):448–454

    Article  Google Scholar 

  • Coleman RG, Gregory RT, Brown GF (1983) Cenozoic volcanic rocks of Saudi Arabia. US Geological Survey, Open-File Report, pp 83–788

  • Condit CD, Connor CB (1996) Recurrence rates of volcanism in basaltic volcanic fields: an example from the Springerville volcanic field, Arizona. Geol Soc Am Bull 108(10):1225–1241

    Article  Google Scholar 

  • Connor CB (1990) Cinder cone clustering in the TransMexican Volcanic Belt: implications for structural and petrologic models. J Geophys Res 95(B12):19395–19405

    Article  Google Scholar 

  • Connor CB, Bebbington M, Marzocchi W (2015) Probabilistic volcanic Hazard assessment. In: The Encyclopedia of Volcanoes, pp 897–910

    Chapter  Google Scholar 

  • Connor CB, Condit CD, Crumpler LS, Aubele JC (1992) Evidence of regional structural controls on vent distribution: Springerville volcanic field, Arizona. J Geophys Res 97(B9):12349–12359

    Article  Google Scholar 

  • Connor CB, Connor LJ, Germa A, Richardson JA, Bebbington MS, Gallant E, Saballos A (2019) How to use kernel density estimation as a diagnostic and forecasting tool for distributed volcanic vents. Stat Volcanol 4(3):1–25

    Article  Google Scholar 

  • Connor CB, Connor LJ, Jaquet O, Wallace L, Kiyosugi K, Chapman N, Sparks S, Goto J (2013) Spatial and temporal distribution of future volcanism in the Chugoku region. Nuclear waste Management Organization of Japan (NUMO), NUMO-TR-13-03

    Google Scholar 

  • Connor CB, Conway FM (2000) Basaltic volcanic fields. In: Sigurdsson H, Houghton B, McNutt SR, Rymer H, Stix J (eds) The encyclopedia of volcanoes, pp 331–343

    Google Scholar 

  • Connor CB, Hill BE (1995) Three nonhomogeneous Poisson models for the probability of basaltic volcanism: application to the Yucca Mountain region, Nevada. J Geophys Res 100(B6):10107–10125

    Article  Google Scholar 

  • Connor LJ, Connor CB, Meliksetian K, Savov I (2012) Probabilistic approach to modeling lava flow inundation: a lava flow hazard assessment for a nuclear facility in Armenia. J Appl Volcanol 1:1–19

    Article  Google Scholar 

  • Connor CB, Stamatakos JA, Ferrill DA, Hill BE, Ofoegbu GI, Conway FM, Sagar B, Trapp J (2000) Geologic factors controlling patterns of small-volume basaltic volcanism: Application to a volcanic hazards assessment at Yucca Mountain, Nevada. J Geophysical Res 105:417–432

  • Downs DT, Champion DE, Muffler P, Christiansen RL, Clynne MA, Calvert AT (2020) Simultaneous middle Pleistocene eruption of three widespread tholeiitic basalts in northern California (USA): insights into crustal magma transport in an actively extending back arc. Geology 48(12):1216–1220

    Article  Google Scholar 

  • Duncan RA, Al-Amri AM (2013) Timing and composition of volcanic activity at Harrat Lunayyir, western Saudi Arabia. J Volcanol Geotherm Res 260:103–116

    Article  Google Scholar 

  • Duncan RA, Kent AJ, Thornber CR, Schlieder TD, Al-Amri AM (2016) Timing and composition of continental volcanism at Harrat Hutaymah, western Saudi Arabia. J Volcanol Geotherm Res 313:1–14

    Article  Google Scholar 

  • Duong T (2005) Bandwidth selectors for multivariate kernel density estimation. Ph.D. thesis, School of Mathematics and Statistics, University of Western Australia

    Google Scholar 

  • El Difrawy MA, Runge MG, Moufti MR, Cronin SJ, Bebbington M (2013) A first hazard analysis of the quaternary Harrat Al-Madinah volcanic field, Saudi Arabia. J Volcanol Geotherm Res 267:39–46

    Article  Google Scholar 

  • Evarts RC, Conrey RM, Fleck RJ, Hagstrum JT (2009) The boring volcanic field of the Portland-Vancouver area, Oregon and Washington: tectonically anomalous forearc volcanism in an urban setting. Field Guides 15:253–270

    Google Scholar 

  • Gallant E, Cole L, Connor C, Donovan A, Molisee D, Morin J, Walshe R, Wetmore P (2021) Modelling eruptive source parameters in distributed volcanic fields. Volcanica 4(2):325–343

    Article  Google Scholar 

  • Gallant E, Richardson J, Connor C, Wetmore P, Connor L (2018) A new approach to probabilistic lava flow hazard assessments, applied to the Idaho National Laboratory, eastern Snake River plain, Idaho, USA. Geology 46(10):895–898

    Article  Google Scholar 

  • Gonnermann H, Taisne B (2015) Magma transport in dikes. The Encyclopedia of Volcanoes 215–224

  • Grosse P, Ramacciotti MLO, Fochi FE, Guzmán S, Orihashi Y, Sumino H (2020) Geomorphology, morphometry, spatial distribution and ages of mafic monogenetic volcanoes of the Peinado and Incahuasi fields, southernmost central volcanic zone of the Andes. J Volcanol Geotherm Res 401:106966

    Article  Google Scholar 

  • Haag MB, Baez WA, Sommer CA, Arnosio JM, Filipovich RE (2019) Geomorphology and spatial distribution of monogenetic volcanoes in the southern Puna plateau (NW Argentina). Geomorphology 342:196–209

    Article  Google Scholar 

  • Hardee HC (1982) Incipient magma chamber formation as a result of repetitive intrusions. Bull Volcanol 45:41–49

    Article  Google Scholar 

  • Heath M, Phillips D, Matchan EL (2020) Basalt lava flows of the intraplate newer Volcanic Province in south-East Australia (Melbourne region): 40Ar/39Ar geochronology reveals~ 8 ma of episodic activity. J Volcanol Geotherm Res 389:106730

    Article  Google Scholar 

  • Ho CH, Smith EI, Feuerbach DL, Naumann TR (1991) Eruptive probability calculation for the Yucca Mountain site, USA: statistical estimation of recurrence rates. Bull Volcanol 54:50–56

    Article  Google Scholar 

  • Hopkins JL, Smid ER, Eccles JD, Hayes JL, Hayward BW, McGee LE, Wijk KV, Wilson TM, Cronin SJ, Leonard GS, Lindsay JM, Németh K, Smith IE (2020) Auckland volcanic field magmatism, volcanism, and hazard: a review. N Z J Geol Geophys 64(2–3):213–234

    Google Scholar 

  • Karlstrom L, Dufek J, Manga M (2009) Organization of volcanic plumbing through magmatic lensing by magma chambers and volcanic loads. J Geophys Res 114(B10204):1–16

  • Kear D (1964) Volcanic alignments north and west of New Zealand's central volcanic region. N Z J Geol Geophys 7:24–44

    Article  Google Scholar 

  • Kereszturi G, Németh K, Moufti MR, Cappello A, Murcia H, Ganci G, Negro CD, Procter J, Zahran HMA (2016) Emplacement conditions of the 1256 AD Al-Madinah lava flow field in Harrat Rahat, Kingdom of Saudi Arabia — insights from surface morphology and lava flow simulations. J Volcanol Geotherm Res 309:14–30

    Article  Google Scholar 

  • Konrad K, Graham DW, Thornber CR, Duncan RA, Kent AJ, Al-Amri AM (2016) Asthenosphere–lithosphere interactions in Western Saudi Arabia: inferences from 3He/4He in xenoliths and lava flows from Harrat Hutaymah. Lithos 248-251:339–352

    Article  Google Scholar 

  • Koppers AA (2002) ArArCALC—software for 40Ar/39Ar age calculations. Comput Geosciences 28(5):605–619

  • Kósik S, Bebbington M, Németh K (2020) Spatio-temporal hazard estimation in the central silicic part of Taupo volcanic zone, New Zealand, based on small to medium volume eruptions. Bull Volcanol 82:1–15

    Article  Google Scholar 

  • Kugaenko Y, Volynets AO (2019) Magmatic plumbing systems of the monogenetic volcanic fields: a case study of Tolbachinsky Dol, Kamchatka. J Volcanol Geotherm Res 383:63–76

    Article  Google Scholar 

  • Kuiper KF, Deino A, Hilgen FJ, Krijgsman W, Renne PR, & Wijbrans JR (2008) Synchronizing rock clocks of Earth history. Sci 320(5875):500–504

  • Kuntz, M.A., Covington, H.R., and Schorr, L.J., 1992. An overview of basaltic volcanism of the eastern Snake River plain, Idaho, in link, P.K., et al., eds., Regional Geology of Eastern Idaho and Western Wyoming: Geological Society of America Memoirs, 179, 227–268

  • Le Corvec N, Spörli KB, Rowland J, Lindsay J (2013) Spatial distribution and alignments of volcanic centers: clues to the formation of monogenetic volcanic fields. Earth Sci Rev 124:96–114

    Article  Google Scholar 

  • Lim JA, Chang SJ, Mai PM, Zahran H (2020) Asthenospheric flow of plume material beneath Arabia inferred from S wave traveltime tomography. Journal of geophysical research: solid. Earth 125(8):e2020JB019668

    Google Scholar 

  • Lindsay JM, Moufti MR (2014) Assessing volcanic risk in Saudi Arabia. Eos Trans Am GeophysUnion 95(31):277–278

    Article  Google Scholar 

  • Lutz TM (1986) An analysis of the orientations of large-scale crustal structures: a statistical approach based on areal distributions of pointlike features. J Geophys Research 91(B1):421–434

    Article  Google Scholar 

  • Lutz TM, Gutmann JT (1995) An improved method for determining and characterizing alignments of pointlike features and its implications for the pinacate volcanic field, Sonora, Mexico. J Geophys Res 100(B9):17659–17670

    Article  Google Scholar 

  • Miyabuchi Y, Hanada D, Niimi H, Kobayashi T (2013) Stratigraphy, grain-size and component characteristics of the 2011 Shinmoedake eruption deposits, Kirishima volcano, Japan. J Volcanol Geotherm Res 258:31–46

    Article  Google Scholar 

  • Mnich ME, Condit CD (2018) Basaltic magmatic mapping: a suggested methodology and the resulting petrologic and volcanic hazard implications from the Springerville volcanic field, east Central Arizona. J Volcanol Geotherm Res 366:58–73

    Article  Google Scholar 

  • Morfulis M, Báez W, Retamoso S, Bardelli L, Filipovich R, Sommer CA (2020) Quantitative spatial distribution analysis of mafic monogenic volcanism in the southern Puna, Argentina: implications for magma production rates and structural control during its ascent. J S Am Earth Sci 104:102852

    Article  Google Scholar 

  • Muffler LJP, Clynne MA, Calvert AT, Champion DE (2011) Diverse, discrete, mantle-derived batches of basalt erupted along a short normal fault zone: the poison Lake chain, southernmost cascades. Bulletin 123(11–12):2177–2200

    Google Scholar 

  • Murcia H, Lindsay JM, Németh K, Smith IEM, Cronin SJ, Moufti MRH, El-Masry NN, Niedermann S (2017) Geology and geochemistry of Late Quaternary volcanism in northern Harrat Rahat, Kingdom of Saudi Arabia: implications for eruption dynamics, regional stratigraphy and magma evolution. Geol Soc Lond, Spec Publ 446:173–204

    Article  Google Scholar 

  • Nieto-Torres A, Martin Del Pozzo AL (2019) Spatio-temporal hazard assessment of a monogenetic volcanic field, near México City. J Volcanol Geotherm Res 371:46–58

    Article  Google Scholar 

  • O’Leary DW, Mankinen EA, Blakely RJ, Langenheim VE, Ponce DA (2002) Aeromagnetic expression of buried basaltic volcanoes near Yucca Mountain, Nevada, USGS-OFR-02-020

    Google Scholar 

  • Pallister JS, McCausland WA, Jónsson S, Lu Z, Zahran HM, El Hadidy S, Aburukbah A, Stewart LCF, Lundgren PR, White RA, Moufti MRH (2010) Broad accommodation of rift-related extension recorded by dyke intrusion in Saudi Arabia. Nat Geosci 3(10):705

    Article  Google Scholar 

  • Richardson JA, Bleacher JE, Glaze LS (2013) The volcanic history of Syria Planum, Mars. J Volcanol Geotherm Res 252:1–13

    Article  Google Scholar 

  • Roobol, M.J., Camp, V.E., 1991. Geology of the Cenozoic lava fields of Harrats Khaybar, Ithnayn and Kura, Kingdom of Saudi Arabia: Saudi Arabian Directorate General of Mineral Resources Geoscience Map GM-131, scale 1:250,000 (in press)

    Google Scholar 

  • Runge MG, Bebbington MS, Cronin SJ, Lindsay JM, Kenedi CL, Moufti MRH (2014) Vents to events: determining an eruption event record from volcanic vent structures for the Harrat Rahat, Saudi Arabia. Bull Volcanol 76(3):804

    Article  Google Scholar 

  • Runge MG, Bebbington MS, Cronin SJ, Lindsay JM, Moufti MR (2016) Integrating geological and geophysical data to improve probabilistic hazard forecasting of Arabian shield volcanism. J Volcanol Geotherm Res 311:41–59

    Article  Google Scholar 

  • Siebert L, Simkin T, Kimberly P (2011) Volcanoes of the world. University of California Press

  • Sieron K, Juárez Cerrillo SF, González-Zuccolotto K, Córdoba-Montiel F, Connor CB, Connor L, Tapia-McClung H (2021) Morphology and distribution of monogenetic volcanoes in Los Tuxtlas volcanic field, Veracruz, Mexico: implications for hazard assessment. Bull Volcanol 83(7):1–17

    Article  Google Scholar 

  • Stelten ME, Downs DT, Champion DE, Dietterich HR, Calvert AT, Sisson TW, Mahood GA, Zahran H (2020) The timing and compositional evolution of volcanism within northern Harrat Rahat, Kingdom of Saudi Arabia. GSA Bull 132(7–8):1381–1403

    Article  Google Scholar 

  • Stern RJ, Johnson P (2010) Continental lithosphere of the Arabian plate: a geologic, petrologic, and geophysical synthesis. Earth Sci Rev 101(1–2):29–67

    Article  Google Scholar 

  • Sulpizio R, Folch A, Costa A, Scaini C, Dellino P (2012) Hazard assessment of far-range volcanic ash dispersal from a violent Strombolian eruption at Somma-Vesuvius volcano, Naples, Italy: implications on civil aviation. Bull Volcanol 74(9):2205–2218

    Article  Google Scholar 

  • Thomson BJ, Lang NP (2016) Volcanic edifice alignment detection software in MATLAB: test data and preliminary results for shield fields on Venus. Comput Geosci 93:1–11

    Article  Google Scholar 

  • Valente F, Cruz JV, Pimentel A, Coutinho R, Andrade C, Nemésio J, Cordeiro S (2022) Evaluating the impact of explosive volcanic eruptions on a groundwater-fed water supply system: an exploratory study in Ponta Delgada, São Miguel (Azores, Portugal). Water 14(7):1022

    Article  Google Scholar 

  • Valentine GA (1993) Note on the distribution of basaltic volcanism associated with large silicic centers. J Volcanol Geotherm Res 56(1–2):167–170

    Article  Google Scholar 

  • Valentine GA, Connor CB (2015) Basaltic Volcanic Fields. In: The Encyclopedia of Volcanoes, pp 423–439

    Chapter  Google Scholar 

  • Valentine GA, Ort MH, Cortés JA (2021) Quaternary basaltic volcanic fields of the American southwest. Geosphere 17(6):2144–2171

    Article  Google Scholar 

  • Von Veh MW, Németh K (2009) An assessment of the alignments of vents based on geostatistical analysis in the Auckland volcanic field, New Zealand. Geomorphologie 15(3):175–186

    Article  Google Scholar 

  • Weller JN (2004) Bayesian inference in forecasting volcanic hazards: an example from Armenia. M.S. Thesis, University of South Florida

  • Wu J, Cronin SJ, Rowe MC, Wolff JA, Barker SJ, Fu B, Boroughs S (2021) Crustal evolution leading to successive rhyolitic supereruptions in the Jemez Mountains volcanic field, New Mexico, USA. Lithos 396:106201

    Article  Google Scholar 

  • Wu J, Rowe MC, Cronin SJ, Wolff JA, Fu B (2022) Long-lived dacitic magmatic systems and recharge dynamics in the Jemez Mountains volcanic field, western USA. Contrib Mineral Petrol 177(6):1–17

    Article  Google Scholar 

  • Zahran HM, Stewart ICF, Johnson PR, Basahel MH (2003) Aeromagnetic-anomaly maps of central and western Saudi Arabia. Saudi Geological Survey Open-File Report SGS-OF-2002-8, 6

  • Zhang D, Lutz T (1989) Structural control of igneous complexes and kimberlites: a new statistical method. Tectonophysics 159(1–2):137–148

    Article  Google Scholar 

Download references


This work constitutes part of the MS thesis of Abdullah Alohali and is supported by the Saudi Arabian Government and Oregon State University. Special thanks to the Saudi Geological Survey (SGS) for providing digital maps of the coalesced harrats (after Camp et al. 1991), Marcos Morfulis for his practical suggestions in analyzing and visualizing alignments. Alohali is extremely grateful to his colleague, Suhail Alhejji, for his continued support and guidance, Eyad Moria for his help in generating random variates and in the probability distribution, and Dr. Elkhedr Ibrahim from King Saud University for providing Harrat Khaybar aeromagnetic data. Cronin thanks King Abdullah Aziz University for support of field work and logistics in Saudi Arabia. Robert Duncan acknowledges support from the US NSF Award# EAR 1524675 for the geochronological analyses presented in this paper and thanks Prof. Abdullah Alamri at King Said University for logistical support. We are also grateful to the associate editor, Charles B. Connor, and the two anonymous reviewers for their constructive comments and thoughtful reviews that substantially improved this manuscript.


This research was gratefully supported by Researchers Supporting Project number (RSP2022R432), King Saud University, Riyadh, Saudi Arabia.

Author information

Authors and Affiliations



Abdullah Alohali: Conceptualization, Methodology, Formal analysis, Investigation, Resources, Writing – Original Draft, Writing – Review & Editing, Visualization, and Funding acquisition. Daniel Bertin: Methodology, Formal analysis, Writing – Review & Editing, and Visualization. Shanaka L. de Silva: Conceptualization, Writing – Review & Editing, Supervision, Project administration, and Funding acquisition. Shane Cronin: Data acquisition, data analysis, and Writing – Review & Editing. Robert Duncan: Data acquisition, data analysis, and Writing – Review & Editing. Saleh Qaysi: Investigation, Resources, Writing – Review & Editing, and Funding acquisition. Mohammed Rashad Hassan Moufti: Field guidance, review and Editing. The author(s) read and approved the final manuscript.

Corresponding author

Correspondence to Abdullah Alohali.

Ethics declarations

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1.

Excel file containing separate sheets for each dataset analyzed in this study. It includes spatial and temporal datasets for Harrat Khaybar such as vent and fissure locations, bandwidth matrices, age data, recurrence rates, and spatio-temporal cumulative probabilities.

Additional file 2.

Analytical data for each of the 15 age determinations of Harrat Khaybar rock that was used in this study.

Additional file 3.

Spatial probability isocontour maps of age group 1, 2, 3 and 5. (a, b, c, and d) Maps depict the distribution of probabilities throughout the harrat based on the locations of visible vents in age group 1, 2, 3 and 5, respectively, without using age-weighting procedure.

Additional file 4.

Reduced-to-pole aeromagnetic maps of Harrat Khaybar. It shows the inconsistency between the visible volcanic vent locations and magnetic anomalies. a) Raw aeromagnetic data. b) Aeromagnetic data clipped to the extent of Harrat Khaybar lavas. c) Stratigraphic unit boundaries of Harrat Khaybar plotted over the aeromagnetic data, corresponding to the boundaries in Fig. 2a. d) Surface vent locations plotted over the aeromagnetic data. The scale in a) represents magnetic intensity.

Additional file 5.

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

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Alohali, A., Bertin, D., de Silva, S. et al. Spatio-temporal forecasting of future volcanism at Harrat Khaybar, Saudi Arabia. J Appl. Volcanol. 11, 12 (2022).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI:


  • Harrat Khaybar
  • Arabian Harrats
  • Eruption event record
  • Spatio-temporal analysis
  • Basaltic volcanic field
  • Distributed volcanism
  • Intraplate volcanism