We present a detailed geologic investigation of Pleistocene to Holocene mafic volcanism within the northernmost part of the Harrat Rahat volcanic field, proximal to the city of Al-Madinah, Saudi Arabia. Our study area covers ∼570 km2, and encompasses lava flows, scoria cones, and shield volcanoes of 32 mapped eruptive units consisting of continental, intraplate alkalic and tholeiitic basalts, hawaiites, and a mugearite that erupted from at least 1014 ± 14 ka to a single Holocene event at 1256 A.D. Typical lava flows are roughly 10–15 km long, although they reach nearly 23 km, 1–3 km wide, and ∼10 m thick. The majority of eruptives in our study area erupted ca. 400–340 ka and ca. 180–100 ka. Despite small individual volumes (<1 km3 dense rock equivalent), each unit resulted from eruption of a distinct magma batch that was influenced by clinopyroxene, olivine, and plagioclase fractionation. Some of these units are interpreted to have undergone magma mixing prior to eruption. Combining our age determinations, geochemistry, and paleomagnetic data sets indicates that several eruptions were temporally and/or spatially clustered. Aligned scoria cones and elongate vent edifices were constructed atop fissure vent systems that reflect the stress field controlling dike ascent through the middle to upper crust.

Small-volume mafic volcanoes (rarely exceeding ∼1 km3) are the most prevalent continental volcanic landform on Earth. These landforms have been documented within volcanic fields associated with subduction zones, rift systems, hot spots, and intraplate settings (e.g., de Silva and Lindsay, 2015; Valentine and Connor, 2015) with >200 volcanic fields considered active during the Holocene (Siebert et al., 2010). Eruptive activity is often spatially and temporally (within hundreds to tens of thousands of years) clustered, with tectonism or preexisting structures controlling or strongly influencing vent locations (e.g., Condit and Connor, 1996; Muffler et al., 2011; Le Corvec et al., 2013; Fleck et al., 2014; Deligne et al., 2016). Spatial, temporal, and geochemical evolution of some small-volume mafic volcanic fields have been studied in considerable detail (e.g., Kuntz et al., 1986; Conway et al., 1997; Shaw et al., 2003; Valentine and Perry, 2007; Duncan and Al-Amri, 2013; Fleck et al., 2014; Deligne et al., 2016; Duncan et al., 2016), but most have received little attention. In addition, few investigations of individual eruptions, and their associated hazards, within volcanic fields have been undertaken (e.g., Valentine and Gregg, 2008; Rowland et al., 2009; Brand et al., 2014). Investigations into the evolution of mafic volcanic fields yield insights on eruptive styles, volumes, and durations, as well as the episodic nature of volcanic activity.

The western Arabian plate encompasses at least 15 continental, intraplate volcanic fields (known in Arabic as harrat) that stretch >3000 km south to north from Yemen through Saudi Arabia to Jordan, Syria, and Turkey (Figs. 1A and 1B). In total, these volcanic fields comprise one of the largest alkalic volcanic provinces on Earth, covering an area of ∼180,000 km2 (Coleman et al., 1983). Historic activity has been recorded with at least 21 eruptions proposed over the past ∼1500 years, the last of which occurred at Dhamar (within the Yemen Traps in Fig. 1B) in northern Yemen in 1937 A.D. (Ambraseys et al., 2005; Siebert et al., 2010). Eruptions have also occurred near the southern end of the Red Sea rift as recently as 2013 with fatalities reported during a 2007 eruption. However, most proposed Holocene eruptions are both unconfirmed and uninvestigated. In addition, many of these volcanic fields remain poorly studied, but growing populations and infrastructure across the region increase the risks from potential future activity (e.g., Lindsay and Moufti, 2014). As a result, the timing, composition, and magmatic processes controlling and influencing activity within these late Cenozoic volcanic fields are being investigated (e.g., Shaw et al., 2003; Moufti et al., 2012, 2013; Duncan and Al-Amri, 2013; Murcia et al., 2015, 2017; Abdelwahed et al., 2016; Duncan et al., 2016; Konrad et al., 2016).

Most of these volcanic fields are situated in sparsely populated regions, but one of the Holocene eruptions occurred within the northernmost part of the Harrat Rahat volcanic field (Fig. 2), near the city of Al-Madinah Al-Munawarah (hereafter referred to as Al-Madinah; Figs. 3A–3D). Only a few individual eruptions from Harrat Rahat have been studied in detail, with particular emphasis on the youngest appearing scoria cones and lava flows within the volcanic field’s more northerly confines close to Al-Madinah (e.g., Camp et al., 1987; Murcia et al., 2015, 2017). All remaining studies have dealt with this volcanic field in its entirety and have primarily focused on magmatic processes (e.g., Camp and Roobol, 1989; Moufti et al., 2012, 2013). Still needed is a systematic study of eruptive strata for individual volcanoes and volcanic centers within Harrat Rahat. This is critical for: (1) understanding the overall spatial, temporal, and compositional evolution of the volcanic field; (2) determining the timescales of magmatic processes in the mantle and crust; and (3) understanding hazards and risks associated with the varied styles of volcanism in the region.

Here we use field relations, petrography, geochemistry, paleomagnetism, and 40Ar/39Ar radiometric and 36Cl cosmogenic surface-exposure dating methods to discriminate individual mafic eruptive units within the northernmost part of Harrat Rahat. These results are used to understand the timing, composition, and eruptive processes of mafic volcanism partly obscured by anthropogenic processes of urbanization.

Harrat Rahat is the largest volcanic field within Saudi Arabia at 50–75 km wide, ∼300 km long (Fig. 2), covering an area of ∼20,000 km2 with a volume of ∼2000 km3, and encompassing >900 observable vents (Coleman et al., 1983; Camp and Roobol, 1989, 1991; Runge et al., 2014). It is a composite of four smaller volcanic fields (Harrat Ar Rukhq, Harrat Turrah, Harrat Bani Abdullah, and Harrat Rashid; Fig. 2) that have coalesced. The northernmost part of Harrat Rahat is sometimes referred to as Harrat Al-Madinah (Moufti, 1985); however, we use Harrat Rahat throughout. In this and other harrats in western Saudi Arabia, dominantly basaltic lava flows were sufficiently clustered and voluminous to obscure underlying Precambrian rocks of the Arabian shield, except for scattered kīpukas along the field’s margins. Harrat Rahat is situated within the western part the Arabia plate in the center of a northward-aligned series of volcanic fields (Figs. 1A and 1B). On a more local scale Harrat Rahat itself has a north to south trend, aligning with Harrat Khaybar, Kura, and Ithnayn to the north (separated by ∼50 km), along a feature termed the Makkah-Madinah-Nafud volcanic line (Camp and Roobol, 1992). Clustering of vents along a north to south lineament in northern Harrat Rahat constructed a topographic crest that stands 100–400 m above the adjacent Precambrian basement, with the result that most of the exposed lavas flowed either east or west from this main vent axis (Fig. 2). The northern end of Harrat Rahat is an exception, with flows that also traveled northeast, north, and northwest from the main vent axis into low-lying areas in part occupied by the modern city of Al-Madinah. Volcanism is proposed to have initiated by at least ca. 10 Ma in Harrat Rahat with the most recent eruption occurring in 1256 A.D. immediately southeast of Al-Madinah.

Camp and Roobol (1989, 1991) mapped the volcanic geology of the entire ∼300 km length of Harrat Rahat, mainly by examining and interpreting aerial images, but widespread sampling allowed them to estimate compositions consisting of ∼62% tholeiitic basalt, ∼21% alkalic basalt, ∼13% hawaiite, and ∼4% more evolved compositions of mugearite, benmoreite, and trachyte (nomenclature after Cox et al., 1979). Tholeiitic basalts dominated early volcanism producing now deeply eroded lava flows, some of which extended west to the Red Sea coastal plain, along with initially minor volumes of alkalic basalts. Alkalic basalts replaced tholeiites as the dominant eruptives during the Pliocene and have been contemporaneous with minor volumes of evolved compositions (i.e., hawaiite, mugearite, benmoreite, and trachyte) since the Pleistocene. Basalt, hawaiite, and mugearite magmas all predominantly vented in Hawaiian and Strombolian eruptions that built scoria cones and lava flows that can reach >20 km length. Maar craters and fine-grained mafic tuff rings indicate phreatomagmatic activity through interaction with groundwater (e.g., Camp and Roobol, 1989; Murcia et al., 2015, 2017), although such features are restricted to the high-standing main vent axis at >10 km south of Al-Madinah. Air-fall tephra deposits are rarely preserved, which is likely the result of their small volumes, coupled with a dry, unvegetated environment where eolian and fluvial processes quickly remobilize fine-grained unconsolidated materials. The only well-preserved and studied air-fall tephra is a near-vent (up to ∼5 km distance) deposit created by fire fountaining during the 1256 A.D. eruption (Kawabata et al., 2015). Benmoreite and trachyte magmas tended to erupt as Peléan domes and spines, including common cryptodomes, before occasionally collapsing and generating small-volume pyroclastic flows consisting dominantly of poorly inflated juvenile clasts. Pumiceous tephra air-fall deposits, as well as pumiceous pyroclastic flows, are rare but not entirely absent (e.g., Moufti, 1985; Camp and Roobol, 1989). Benmoreite and trachyte eruptives are only present >10 km south of our study area and as such will not be discussed further.

Prior to this study, late Cenozoic volcanic stratigraphy of Harrat Rahat was divided into the Shawahit (ca. 10–2.5 Ma), Hammah (ca. 2.5–1.7 Ma), and Madinah basalts (ca. 1.7 Ma to present) by Camp and Roobol (1989, 1991). The Madinah basalts were further subdivided into seven Quaternary units (Qm1–7) based on geomorphic criteria, such as the degree of erosion and size and abundance of loess-filled depressions. Absolute ages within the Qm1–7 system were based on eight K-Ar radiometric dates from northern Harrat Rahat, as well as the historical account of the 1256 A.D. eruption (Camp and Roobol, 1989). These age categories have been supplanted by 25 40Ar/39Ar dates from Moufti et al. (2013). Based on their 40Ar/39Ar ages, Moufti et al. (2013) redefined the Qm1–7 system as lasting from 10.31 Ma to the present. However, each of the individual subdivisions overlap using the new stratigraphic age constraints suggested by Moufti et al. (2013). In addition, those K-Ar and 40Ar/39Ar age determinations are commonly contradicted by stratigraphic superposition. Only one of the previous ages is from a lava flow within our study area. That 40Ar/39Ar age is reported in Moufti et al. (2013) as 270 ± 60 ka for the unit we refer to as the basalt of al Anahi 3 (ban3: unit abbreviations in bold and defined in Table 1). Thus, a lack of age determinations on mafic eruptives within and close to Al-Madinah and inconsistencies between published radiometric ages and stratigraphic superposition highlight a need for reliable age and geologic constraints throughout the northernmost part of Harrat Rahat, and within Harrat Rahat in general. We present 19 new age determinations on 32 mapped eruptive units in a ∼570 km2 area in the northwesternmost part of the Harrat Rahat volcanic field, which includes the city of Al-Madinah (Figs. 4 and 5). Many of the contacts depicted by Camp and Roobol (1991) have been identified here as internal lava-flow structures (i.e., levees, channels, etc.).

The most prominent volcanic landforms are the rough surfaces of lava flows and their source vents, which include scoria cones and shield volcanoes (Figs. 3A–3D). Where possible, vents have been correlated to their associated lava flows, but some vents remain unidentified due to scoria cone collapse and burial by younger lava flows. Most cones are steep sided and rise >100 m above their surrounding lava flows, although some of the smallest volume eruptions have scoria cones only a few tens of meters high. Many eruptions created scoria cones with nested craters or clusters of north- to northwest-aligned craters or cones. The latter result from fissure eruptions, such as the ∼2.25-km-long fissure that formed seven north-northwest–aligned scoria cones during the 1256 A.D. basalt of al Labah (bla) eruption (Figs. 3A and 4). Other examples of elongate scoria cone complexes interpreted as fissure eruptions are evident on Figure 4, such as the basalt of al Du’aythah (bdu), basalt of ad Duwaykhilah (bdw), basalt of as Suddiyah (bsu), and basalt of Umm Nathilah (bun). These scoria cones also typically display breaches where lava drained out, and rafted cone material was carried along by the lava flows, sometimes to distances >1 km from vent. In general, vents are composed of vesicular, vitric bombs (occasional spindle bombs) and lapilli, as well as spatter and agglutinate. Due to loose consolidation, many of the older (>100 k.y.) scoria cones are heavily eroded, exposing deeply oxidized, syneruptive reddish-brown material. Only the youngest scoria cones retain a black, vitric uneroded appearance. Lava flows in northernmost Harrat Rahat can attain thicknesses of >10 m and produce near-vent pāhoehoe, broad ‘a‘ā channels with levees, and transitional morphologies (Figs. 3A, 3B, and 3D). Small shield volcanoes are sparse within northern Harrat Rahat, and only two are present near Al-Madinah—the vents for the basalt of al Huzaym (bhu) and basalt of the Reverse (bre; Fig. 4). These shields have gentle slopes of only a few degrees and can measure up to ∼4 km in diameter with interior craters of ∼1 km in diameter.

The majority of eruptives near Al-Madinah are alkalic basalts with minor volumes of tholeiitic basalts, hawaiites, and a single unit of mugearite composition (Table 1; Fig. 6A). In addition, there is both petrographic and geochemical evidence for mixing of distinctive basaltic components prior to a few eruptions (e.g., bla in Camp et al., 1987). All Harrat Rahat mafic eruptives have broadly similar petrographic characteristics and mineral assemblages from aphyric to olivine (± chromian spinel inclusions) ± plagioclase phenocrysts; and olivine + plagioclase + pyroxene + Fe-Ti oxide + glass-bearing groundmass. Normally zoned olivine phenocrysts range from ∼1%–10%, 1–5 mm in size (with rare phenocrysts to >20 mm), and have a variety of colors (dark to pale green, honey yellow, and varying degrees of iddingsite) and textures from resorbed to partially resorbed, corroded, and skeletal (Figs. 7A–7D). Plagioclase phenocrysts (typically labradorite with minor bytownite) show a similar range of characteristics, ranging from <1%–5% (can achieve >50 mm in size but are mostly <10 mm), and many are euhedral to resorbed with overgrowth rims (Figs. 7A–7D). Pyroxene phenocrysts are rare within mafic lavas, but instead clinopyroxene is almost exclusively confined to interstitial spaces within groundmass as augite, diopside, and hedenbergite. Equant to subequant magnetite and titanomagnetite (from 10 to >100 µm) are ubiquitous throughout the groundmass of Harrat Rahat eruptives, and the more evolved compositions (low-MgO alkalic basalt, hawaiite, and mugearite) contain trace amounts of apatite and amphibole (i.e., kaersutite, ferro-kaersutite). At least eight of the basaltic eruptives (bla, bdw, bre, bhu, brm—basalt of Red Mountain; bra—basalt of al Rafi’ah; bqr—basalt of Quraydah; and bmu—basalt of al Mustarah) contain dark greenish-brown chromian spinels (i.e., picotite) hosted in olivine phenocrysts (Figs. 7B–7D), of which some have magnetite or titanomagnetite rims that may have grown during brief stalling or slow ascent rates (i.e., Kamenetsky et al., 2001).

Shield volcanoes are typically tholeiitic with a coarse-grained, diktytaxitic groundmass, although some are holocrystalline with seriate textures (Fig. 7A). Tholeiitic basalts range from aphyric to phenocryst rich. Olivine is typically the only phenocryst phase at ∼5%–10%, with skeletal, corroded, or euhedral textures, and colors ranging from shades of green, honey yellow, or heavily iddingsitized. Groundmass plagioclase laths are ubiquitous, but plagioclase phenocrysts are rare in Harrat Rahat tholeiitic basalts.

Physical characteristics for individual units vary. Table 1 summarizes estimates for each mapped unit’s estimated area, volume, and flow length from vent. Due to the severity of anthropogenic modification of vents and lava flows around Al-Madinah, erupted volumes were calculated using a uniform lava-flow thickness of 10 m (although some may get to 15 m thick in areas). This is consistent with field observations in the less disturbed parts of Harrat Rahat further from Al-Madinah. Older lava flows are exposed discontinuously due to partial coverage by younger lava flows, and their limits have been conservatively reconstructed where possible. Minimum estimates of scoria cone volumes were made by drawing multiple cross sections to assess average cone heights and multiplying by cone area. These cones are typically strongly elongate in the northwest direction, and many have undergone breaches or syneruptive to posteruptive collapse, precluding the use of a simple geometric formula for the volume of a cone. In addition, with the exception of unit bdu, scoria cones comprise ≤1% of the total dense rock equivalent (DRE) volume of eruptions where exposures provide reliable constraints on vent and lava-flow area and volume estimates. Lava-flow DRE volumes are calculated using a lava density equal to that of a basaltic magma density of 2500 kg/m3 (Crosweller et al., 2012) and a scoria cone density of 1500 kg/m3 (McGetchin et al., 1974).

Individual eruptive products range in area from the four small scoria cones of bdu at ∼0.1 km2 to ban3 at ∼69 km2. As expected, these also have the smallest and largest DRE volumes at ∼0.002 and ∼0.7 km3, respectively. Overall, mapped eruptive units cover an average area of ∼21 km2 and extruded an average DRE volume of ∼0.2 km3. However, since older eruptives are buried beneath younger lava flows to the point that reconstructing their total areas and volumes is precluded, we select only the youngest eruptives to provide a more precise estimate of the total area and erupted volumes of individual units, which are relevant to potential future eruptions. Using the 13 youngest eruptions (<200 k.y.) that have vents tied to lava flows gives a higher average area of ∼34 km2 and average DRE volume of ∼0.3 km3. Lava-flow lengths vary from slightly more than 20 km, with bla reaching ∼23 km in length, to the short, stubby lava flows of bdu that are only 115 and 200 m long. The same 13 eruptions yield an average lava-flow length of ∼13.5 km. This only takes into account the length of the longest lobe of lava from the vent, despite some units containing multiple flows extending down multiple drainage channels. For example, ban3 has a lobe of lava that reaches 20.6 km northwest of its vent, but also one that reaches 14.3 km west-northwest, and a third that reaches 15 km west-southwest. These more southerly two lava flows are not shown on Figure 4 as they are beyond the limits of the map area.

Discrimination of Map Units

We mapped individual eruptive units in the field through outcrop-scale and hand-specimen petrography and observing lava-flow contacts and superposition, supported by 0.5 m resolution light detection and ranging (lidar) and satellite photogrammetric digital topographic bases. Petrographic differences are noticeable between almost all immediately adjacent units, such as the abundance, size, and texture of both phenocrysts and groundmass (Table 1). Thin section petrography, major-oxide and trace-element bulk-rock geochemistry via X-ray fluorescence (XRF) and inductively coupled plasma–mass spectrometry (ICP-MS), paleomagnetic studies, and age determinations via 40Ar/39Ar radiometric and 36Cl cosmogenic surface-exposure dating were obtained later to confirm and further refine field observations. Interrogation of all the aforementioned information obtained through modern geologic mapping techniques identified 32 distinct mafic volcanic units ranging from ca. 1 Ma to 1256 A.D. (Figs. 4 and 5).

No consistent stratigraphic naming scheme has been used within Harrat Rahat; therefore, we assigned new unit names based on landforms documented in topographic maps. Where possible, scoria cone names were used, but for units lacking a scoria cone, names were assigned based on distinguishing geomorphic features or the city district where they were predominantly mapped. Previous assigned names were kept where possible, such as bdu (basalt of al Du’aythah in Murcia et al., 2015). Some previously named scoria cones (i.e., al Anahi) of similar composition are clustered, in which case, we used the previously assigned name and added numeric identifiers with one being the geochronologically oldest (i.e., ban1, ban2, and ban3 for the three scoria cones near al Anahi).

Geochemical Interpretation

Major-oxide and trace-element analyses were performed on 151 mafic lava and scoria samples by XRF, and trace-element analysis on 150 of these samples by ICP-MS at the GeoAnalytical Laboratory at Washington State University in Pullman, Washington. Results of representative samples are presented in Table 2 (see the Supplemental File Table S11 for a complete table of XRF and ICP-MS data and details on the XRF and ICP-MS methods). These data are consistent in major-oxide and trace-element abundances with previously published Harrat Rahat mafic eruptives (Camp et al., 1987; Camp and Roobol, 1989; Murcia et al., 2015, 2017).

Slight to significant differences in geochemistry and petrography (Tables 1 and 2) among immediately adjacent eruptives indicate that each distinct unit represents a discrete magmatic batch or, occasionally, the mixing of two magma batches (e.g., bla, bdw, bhu, and bau—basalt of an Nughayr). Most of the Quaternary volcanic rocks within and around Al-Madinah are alkalic basalts, with minor lava flows and scoria cones of tholeiitic basalt, hawaiite, and mugearite (Fig. 6A). Basaltic eruptives (≥6 wt% MgO) studied here have compositions within the intraplate basalt field of Pearce and Norry (1979) and continental basaltic field of Cabanis and Lecolle (1989), although a few samples plot in the alkalic basalts from intercontinental rift field of the latter (Figs. 6B and 6C). Although Harrat Rahat basalts are generated from partial melting of the mantle, even the most primitive basaltic eruptives are not considered true primary magmas produced solely through partial melting of peridotite (e.g., Camp and Roobol, 1989; Duncan et al., 2016, Konrad et al., 2016). Unmodified primary partial melts in equilibrium with mantle peridotite are distinguished by bulk-rock Mg# values of 65–75, Ni values at >300 ppm, and olivine compositions of at least Fo89 (Roeder and Emslie, 1970; O’Reilly and Griffin, 1984). The most primitive basalts from our analyses contain Mg# values from 58.8 to 66.1 (average of 61.8) and Ni values from 129 to 313 ppm (average of 192 ppm) and have olivine phenocryst compositions that average Fo71. Only a single sample from our data set (R16DD198 from bra) might potentially be considered an unmodified partial melt of peridotite.

Bulk-rock variations are used to distinguish between individual eruptive units, with Figure 8 displaying trends for selected major-oxide, trace-element, and Pearce ratio plots (Pearce, 1968). Mapped units range from 45.17 to 49.22 wt% SiO2 and 3.66–11.47 wt% MgO (Mg# values from 35.3 to 66.1), and 120 of the 151 samples analyzed have ≥6 wt% MgO. Tholeiitic basalts are the most primitive with MgO values of 8.48–11.47 wt% (Mg# values from 58.8 to 66.1). Incompatible components such as TiO2, K2O, and P2O5 are particularly illustrative for correlating units because none of these components are retained within early crystallizing phenocryst phases, and their abundances vary considerably over the limited range of SiO2. For example, K2O varies by a factor of >6 from 0.28 to 1.80 wt%, P2O5 varies by a factor of ∼13 from 0.14 to 1.82 wt%, and TiO2 varies by a factor of ∼3 from 1.26 to 3.68 wt%.

In general, individual eruptive units create plotted arrays indicating concurrent clinopyroxene, olivine, and plagioclase fractionation, as well as crystallization of trace phases (chromian spinel, Fe-Ti oxides, and/or apatite depending on the degree of magma evolution); however, the compositions from each unit typically cluster (Fig. 8). The number of chemical analyses (151) coupled with number of mapped units (32) yield unavoidable compositional overlap, particularly for units interpreted as having undergone magma mixing prior to eruption. Nevertheless, combining geochemistry with petrographic descriptions, field mapping, and stratigraphic superposition for adjacent units supports the accuracy of our mapped contacts. In particular, plotting K2O/TiO2 and K2O/P2O5 versus MgO or SiO2 yields groupings that aid in distinguishing eruptive units, which are analogous to the P/K versus Ti/K Pearce element ratio plot (Fig. 8). These elements are not perfectly incompatible over the compositional range discussed here. However, other incompatible-element pairs (e.g., Zr/Nb in Fig. 8) still support our mapped units not being cogenetic, despite some units having experienced preeruptive or syneruptive magma mixing (e.g., bla, bdw, bhu, and bau).

While developing a petrogenetic model is not the focus of this paper, we can make a few observations. Geochemical variations are strongly influenced by the fractionation, crystallization, and accumulation of phenocryst and trace phases. For example, decreasing Sc and CaO/Al2O3 coincident with decreasing MgO indicate crystal fractionation of Mg-rich clinopyroxene (i.e., David et al., 2000; Smith et al., 2008). The general absence of pyroxene phenocrysts within Harrat Rahat mafic eruptives indicates that this pyroxene grew and separated at depth, perhaps near the base of the crust. Similar scenarios of clinopyroxene fractionation have been proposed in other mafic volcanic fields, including western Arabian plate volcanic fields (e.g., Shaw et al., 2003; Smith et al., 2008; McGee et al., 2013; Murcia et al., 2017). Decreasing SiO2 abundances coincide with decreasing MgO, until MgO falls to ∼6.5 wt%, in agreement with clinopyroxene-dominated fractionation ceasing below ∼6.5 wt% MgO, at which point SiO2 starts to increase (Fig. 8) as a result of olivine and plagioclase becoming the primary fractionating phases due to either a reduction in pressure or the appearance of an amphibole phase. Olivine (± chromian spinel inclusions) accumulation is indicated in some volcanic rocks by high abundances of Ni (129–313 ppm) and Cr (144–640 ppm: Hart and Davis, 1978; Hart and Dunn, 1993), while more evolved compositions show variable accumulation of plagioclase by enrichment of Ba and Sr. Chromian-spinel inclusions within olivine phenocrysts are limited to high-MgO alkalic and tholeiitic basalts. Many chromian spinels have Fe-Ti oxide rims or have been completely transformed to Fe-Ti oxides (particularly when situated near olivine rims or fractures in host phenocrysts). Fractionation of Fe-Ti oxides also begins after MgO falls below ∼6.5 wt%, as TiO2 decreases after this point. Apatite starts to crystallize in some low-MgO alkalic basalts and is present in most hawaiite and mugearite compositions, below ∼6.5 wt% MgO as P2O5 sharply increases.

Samples from 16 units have textures appropriate for 40Ar/39Ar radiometric dating, and samples from three units were collected for 36Cl cosmogenic surface-exposure dating. Details for these experiments are outlined in Tables 3 and 4, and details for the 40Ar/39Ar age determinations are presented in Figure S1 (see the Supplemental File [footnote 1] for 40Ar/39Ar radiometric and 36Cl cosmogenic surface-exposure dating methods and complete 40Ar/39Ar age profiles). New 40Ar/39Ar eruption ages presented here all conform to the observed stratigraphic position of units as mapped in the field. These ages are described below from the oldest in the Early and Middle Pleistocene to the youngest Holocene eruption (as designated by the International Commission on Stratigraphy: http://www.stratigraphy.org/). Tertiary basalts (Tb) that erupted ca. 11 Ma (Camp and Roobol, 1991) cap some of the Precambrian mountains surrounding Al-Madinah as erosional remnants (Fig. 4); however, little is known about their origins and context as related to Harrat Rahat. Therefore, they have not been re-dated as part of this study and are not discussed here.

Early and Middle Pleistocene

The history of northernmost Harrat Rahat begins in the Early Pleistocene. The oldest dated Quaternary volcanic unit within our study area is a small shield volcano that produced bre at 1014 ± 14 ka (all uncertainties are at 1σ unless otherwise stated), located ∼20 km southeast of the city center. This is the only directly dated Early Pleistocene eruptive within our study area (Fig. 9A). The next volcanic unit chronologically, and the oldest radiometrically dated within Al-Madinah itself, is the basalt of the Palms (bpa) at 537.1 ± 9.4 ka. This unit was mapped within the central part of the city, overlain and flanked by younger lava flows. These younger lavas lie to its east and west, forming a shallow valley that contains the center and oldest part of Al-Madinah.

These younger lava flows bordering this valley erupted over a period of ∼90–100 k.y., of which five new radiometric dates provide age constraints. The oldest of this group is the basalt of Rahat (brh) at 411.5 ± 5.9 ka, which is present only as a limited exposure in the southwestern part of the city. This was followed by the basalt of al Urayd (bur) at 395.1 ± 6.2 ka, exposed in the northeastern part of the city. Next to erupt were the basalt of al Harrah al Gharbiyah (bhg) at 376.1 ± 5.7 ka and the basalt of the Powerline Road (bpr) at 372.5 ± 9.8 ka. The former constitutes the western margin for the downtown area, while the latter is a kīpuka of small areal extent (∼0.4 km2) completely surrounded by the 1256 A.D. bla lava flow. The youngest dated unit within this group is bmu at 339.0 ± 4.7 ka, which is a lava flow that defines the northeastern margin of the downtown area. Several other mapped units (bqr and bis—basalt of al Iskan and bqu—basalt of Quba) may fall into this age group based on their stratigraphic superposition between age-constrained eruptives, but heavy urbanization precluded obtaining samples suitable for direct dating. The basalt of Hathm (bha) potentially falls within this group as well; however, it could be much older because it underlies bur at 395.1 ± 6.2 ka and sits on no observable underlying units due to its limited exposed area (∼0.7 km2). Therefore, bha could be considerably older than surrounding units, but abundant carbonate- and zeolite-filled vesicles preclude radiometric dating.

Only two units have radiometric ages between ca. 300 and 200 ka. The first is the low-lying shield volcano and associated lava flows that comprise bhu at 280.4 ± 6.5 ka. Based on a general observation of erosional rills on tholeiitic shields, Camp and Roobol (1991) previously interpreted this shield as the oldest volcanic unit beneath Al-Madinah, and they assigned it to the Hammah basalts, which were thought to have erupted between ca. 2.5–1.7 Ma. However, it observably overlies bhg at 376.1 ± 5.7 ka and underlies the basalt of Shuran (bsh) at 188.3 ± 4.2 ka, which supports its measured age. Unit bhu was followed by eruption of the basalt of Banthane Dam (bbd) at 223.4 ± 8.6 ka, which is exposed near the south-central part of Al-Madinah. This unit provides a minimum age constraint on the underlying, undated bqr, as well as a maximum age constraint on all remaining overlying Middle and Late Pleistocene eruptives within the south-central part of Al-Madinah.

The end of the Middle Pleistocene was a period of intense volcanic activity with at least nine eruptions that sent lava flows close to present-day Al-Madinah, of which five have been dated as part of this study. Unit bsh erupted at 188.3 ± 4.2 ka (weighted mean of three groundmass ages; Table 3), which provides a maximum age constraint on the overlying, undated basalt of al Jassah (bjs). Along the eastern to central part of the map area are three basalts with overlapping age determinations. These are the basalt of al Anahi 1 (ban1) at 136.7 ± 5.3 ka, brm at 134.4 ± 7.5 ka, and basalt of Nubala (bnu) at 129.0 ± 8.2 ka. Despite overlapping age uncertainties, they fall within the correct stratigraphic order from oldest to youngest. In addition, brm (134.4 ± 7.5 ka) provides minimum age constraints on three underlying, undated units, which also share an indirect maximum age constraint by overlying bbd (223.4 ± 8.6 ka). These are the basalt of al Muzayyin (bmz), basalt of ar Rummanah (bru), and bra. A few of these units may be Late Pleistocene based on their age uncertainties.

Late Pleistocene

The oldest unequivocal product of a Late Pleistocene eruption is bun at 103.3 ± 8.5 ka, which erupted near the southwestern limits of Al-Madinah. This was followed by an eruption near the eastern city limits of bsu at 84.9 ± 5.7 ka (weighted mean of two groundmass ages; Table 3).

The 36Cl cosmogenic surface-exposure dating technique yielded an age of 51.0 ± 9.0 ka on the young appearing ban3 lava (Table 4). This is significantly younger than an age of 270 ± 60 ka presented by Moufti et al. (2013) from a 40Ar/39Ar experiment on a sample from within our mapped limits of ban3. However, ban3 has a young-appearing morphology, and our age concurs with its stratigraphic position overlying bnu (129.0 ± 8.2 ka) and other similar-aged basaltic lavas inferred from stratigraphic superposition. This result is further supported by an additional constraint of 52.6 ± 4.1 ka (using a 3He cosmogenic surface-exposure date from Murcia et al., 2017) of a young lava flow directly underlying ban3. The 51.0 ± 9.0 ka age from ban3 also provides a minimum age constraint on several underlying and undated, but young-appearing units. These include the mugearite of Hill 1125 (mh11), basalt of al Anahi 2 (ban2), basalt of Pumyra (bpu), and basalt of Sha’ib al Khakh (bsk), which have maximum age constraints from other dated Middle and Late Pleistocene units. We infer that these undated units all erupted <100 ka based on their young appearance, morphology, and stratigraphic positions.

Younger still, bdw erupted at 25.2 ± 15.9 ka from a vent near the eastern limits of Al-Madinah. This basalt overlies unit bsu (84.9 ± 5.7 ka) and is overlain by the historic 1256 A.D. bla eruptives. The large age uncertainty allows for a possible Holocene eruption; however, achieving a more precise age on such a young, low-K2O sample would be difficult by current 40Ar/39Ar methods.

During the end of the Late Pleistocene, four north-northwest–aligned scoria cones (unit bdu), of which the northern two contain short ∼115–200-m-long lava flows, erupted within a western suburb of Al-Madinah >7 km away from any other Harrat Rahat volcanic rocks (Murcia et al., 2015). There has been some confusion over the age of bdu in the literature. An eruption was suggested for either the year 640 or 641 A.D. for bdu based on a historical account that could describe an eruption in western Saudi Arabia around that time. A brief description in Wafa Al-Wafa Bi’Akhbar Dar Al-Mustafa by Al-Samhoody (1486, reprinted and translated to English in 1955) states that during the reign of the second Rashidi Caliphate of Omar Bin al-Khattab, a small fire rose from Harrat Al-Madinah and soon died out. This prompted Camp and Roobol (1989) to propose that scoria cones and short lava flows now designated as unit bdu were products of the eruption that inspired this historical account. However, these same historical accounts seem to present conflicting locations. Accounts indicate that an earthquake destroyed homes in Al-Madinah in 641 A.D. (Ambraseys et al., 2005), and as a result, an eruption was linked to this event. Yet, Ambraseys et al. (2005) also suggested that this same earthquake may have been related to either one or two eruptions that occurred in 640 A.D. within Harrat Uwayrid, which is at least 300 km north of Al-Madinah (Fig. 1B). Adding to the confusion is that these accounts were written centuries after the possible eruption. Therefore, we collected a sample of bdu for 36Cl cosmogenic surface-exposure dating, which yielded an age of 13.3 ± 1.9 ka (Table 4). This unit could have erupted during the Late Pleistocene or early Holocene given its analytical uncertainty.


The only confirmed Holocene eruption in Harrat Rahat is on the southeastern to eastern margin of Al-Madinah. This is the 1256 A.D. bla eruption, which is also the best studied unit within the volcanic field (e.g., Camp et al., 1987; Kawabata et al., 2015). This unit was dated using the 36Cl cosmogenic surface-exposure technique as a test of the method’s robustness, and an age of 0.81 ± 0.24 ka was calculated (Table 4), which compares well with the recorded age of 0.76 ka.

The bla eruption initiated, after five days of progressively stronger earthquakes felt by the residents of Al-Madinah, on Monday, 26 June 1256 A.D. (1 Jumadi Al-Thani 654 A.H. by the Arabic calendar) when a basaltic fissure began erupting ∼20 km southeast of the city center. In total, this eruption lasted 52 days, ending on Sunday, 20 August 1256 A.D. (27 Rajab 654 A.H. by the Arabic calendar). It formed seven scoria cones along a ∼2.25-km-long fissure vent, which propelled 0.0077 km3 of air-fall tephra to >5 km from vent (Kawabata et al., 2015) and extruded 0.45 km3 of lava that flowed ∼23 km from vent to within ∼8 km of the city center. Historical accounts have been recorded in Wafa Al-Wafa Bi’Akhbar Dar Al-Mustafa by Al-Samhoody (1486, reprinted and translated to English in 1955) and Burton (1893). As the best studied eruption and the only one with detailed historical accounts, bla remains an ideal case study for how future eruptive activity within Harrat Rahat will progress.

Another historic eruption (not including the previously discussed 640 or 641 A.D. eruption) has been proposed at 1293 A.D. and continues to be propagated throughout the literature (e.g., Ambraseys et al., 2005; Moufti et al., 2013; Murcia et al., 2015, 2017), although no products have been identified. Both direct physical evidence in the form of a young-looking scoria cone and lava flow and robust historical accounts are lacking for an eruption at 1293 A.D. within Harrat Rahat. In comparison, the eruption of bla in 1256 A.D. is present as a distinctively young, vitric, black scoria cone complex, air-fall tephra apron, and lava flow with well-preserved textures—features corroborated by well-documented historical records. An eruptive event only 37 years after the 1256 A.D. eruption would be just as obvious in the field and in aerial and satellite imagery and documented by our geochronology of all young-appearing lava flows., particularly since young volcanism was concentrated near the northernmost end of Harrat Rahat, proximal to Al-Madinah, which was already a major population center at the time. We find no evidence of an eruption occurring within Harrat Rahat in 1293 A.D. but cannot rule out an eruption in the sparsely populated regions of Harrat Khaybar or Ithnayn, where young-looking basaltic lava flows and air-fall tephra deposits are present, and which could have been distantly visible from Al-Madinah.

We measured paleomagnetic secular variations from mafic lava flows and scoria cones within northernmost Harrat Rahat for the purposes of correlating individual units (Table 5) independent of field, petrographic, and geochemical interpretations. Additionally, paleomagnetic analyses help unravel discrete temporal groupings that could not otherwise be divided based on the precision of modern dating techniques.

Paleomagnetic samples were collected, processed, and interpreted using the methods of McElhinny (1973). We collected eight, 10-cm-long cores at each site in the field using a handheld, gasoline-powered, 2.5 cm coring drill, and oriented exclusively by sun compass. A 2.5-cm-long specimen from each core was measured using an automated cryogenic magnetometer and then was subjected to alternating field (AF) demagnetization to remove secondary components of magnetization. An isothermal component from nearby lightning strikes was the most frequent source of secondary magnetization in these young volcanic rocks. The characteristic direction of remanent magnetization for each site was calculated using Fisher statistics on data from: (1) line fits of data on vector component diagrams; (2) plane fits on equal area diagrams; and (3) mixtures of lines and planes data. We present site mean directions of magnetization for each group in Table 5 with statistics, including 95% confidence limits (α95).

Obtaining accurate paleomagnetic samples from the youngest (<100 k.y.) Harrat Rahat mafic lava flows is difficult because their surfaces tend to be rubbly and uneroded. However, older units have undergone enough dissection that in situ blocks and deeply incised drainage channels into lava-flow interiors could be sampled easily. Drilling in road cuts and drainage channels where erosion had penetrated several meters to the interior of a lava flow was ideal, but eroded surfaces composed of large blocks regarded as having undergone minimal movement since cooling below their Curie temperature were successful as well. Previous experience indicates specific environments where these interior blocks are most likely to have rotated postmagnetization, such as along lava-flow margins, and these were avoided during drilling. Yet, it is impossible (a priori) to preclude minor rotations, and accordingly, dispersion of directions is observed for some groups. This is particularly apparent in the youngest, rubbly lava flows, and for scoria cones, wherein agglutinated spatter was generally the only candidate for drilling. Potential signs of postcooling rotation include larger α95 values after AF treatment and scattered inclination and declination pairs (i.e., mh11, bqr, and brh).

Paleomagnetic Interpretation

We used paleomagnetic data for correlation purposes, with the assumption that each eruptive unit records a single, consistent direction of remanent magnetization acquired during a geologically instantaneous period of time (i.e., days, months, or years). This assumption is based on consistencies in petrographic and geochemical characteristics combined with mapped stratigraphic superposition of individual units and documentation that bla erupted over 52 days. Similar combination of geochemical, petrographic, and paleomagnetic studies from other volcanic fields has yielded results consistent with brief eruptions (e.g., Muffler et al., 2011; Fleck et al., 2014; Deligne et al., 2016). Therefore, we infer that none of these eruptions occurred over protracted time periods (i.e., multiple years to decades).

Most paleomagnetic results from individual units exhibit a single, uniform direction of magnetization (Fig. 10). For example, the 1256 A.D. bla eruption products provide an ideal paleomagnetic case study. Three sites were drilled for analysis, and all yielded inclination declination pairs that overlap within α95 (Fig. 10). Similar tight groupings are present for most units in Table 5 with the average for each unit displayed in bold. In addition, a handful of these units are interpreted as having undergone magma mixing prior to eruption, and as a result their bulk-rock compositions are appreciably heterogeneous, hindering simple chemical correlations between exposures. For example, bhu (280.4 ± 6.5 ka) was previously assigned by Camp and Roobol (1991) to two separate age categories. They assigned the shield volcano to the Hammah basalts between ca. 2.5 and 1.7 Ma, and the lava flow that comprises the northern part of the unit to age category Qm2, which meant it was thought to have erupted between ca. 1.2 Ma and 900 ka. The shield volcano and lava flow differ in chemical composition and texture, including bulk-rock compositions of tholeiitic and alkalic basalts, respectively (Table 2). Such significant differences would normally hinder correlation, particularly in an urban environment where constructional morphology is largely destroyed and exposures are discontinuous. However, paleomagnetic secular variation from two sites encompassed by the shield volcano and three sites within the lava flows yield the same unique average inclination of 45° and declination of 24° (Table 5; Fig. 11A). Therefore, it is statistically unlikely that these represent two temporally separate eruptions. Moreover, recent demolition of the cone in the middle of the shield volcano has uncovered alkalic basalt similar in composition to the lava flows north of the shield volcano. Thus, paleomagnetic data demonstrate this to have been a single, geologically instantaneous eruption in which an alkalic basalt erupted with lavas flowing north, followed by eruption of a tholeiitic basalt from the same vent, creating the shield volcano, with mixing between disparate magmas prior to and during eruption.

The paleomagnetic results also help to constrain estimates of eruption ages. Several units have paleomagnetic secular variations considered excursional, which is taken here as any virtual geomagnetic pole with a co-latitude >23°. The oldest of these is unit bre (1014 ± 14 ka), which is both reverse polarity and excursional. This unit may have been emplaced during the termination of the Jaramillo Subchron at ca. 1 Ma to 990 ka (Fig. 9A). Similarly, a single scoria cone (Quaternary vents [unit Qv] site R16DC003 in Table 5) geographically located between bbd and bpa (Fig. 4) is also reverse polarity, indicating that it likely erupted sometime between ca. 990 and 780 ka during the Matuyama Chron (Fig. 9A). These are the only reverse polarity units identified within our study area, and they are also interpreted as the oldest. While normal and reverse time intervals can last for >100 k.y., excursional time intervals typically exist over short durations of ≤2 k.y. and coincide with increased rates of geomagnetic secular variation (Champion and Shoemaker, 1977). Several units, apart from bre, have larger values of paleomagnetic secular variation indicating that they erupted during excursional time intervals. Two sites drilled within brh (411.5 ± 5.9 ka) yield inclinations of –13.8° and 1.1° and declinations of 56.7° and 351.4° (Table 5), which is in the right time interval to coincide with an unnamed mid-Brunhes Chron excursion (Oda, 2005).

Some mapped eruptives, while separated geographically and chemically distinct, are temporally clustered. Units bhg (376.1 ± 5.7 ka) and bpr (372.5 ± 9.8 ka) have overlapping age uncertainties at 1σ, suggesting they have a shared eruptive time interval on the order of no more than 10–20 k.y. Remanent magnetic mean directions for each of these units also overlap and are separated by only 4°, with bhg yielding an average inclination of 55° and declination of 347° and bpr having an average inclination of 54° and declination of 354° (Fig. 11A). Based on usual rates of geomagnetic secular variation (Champion and Shoemaker, 1977; Merrill and McFadden, 1999), these eruptions are interpreted to have been separated by at most a few centuries. Similarly, the aforementioned unit of bhu (280.4 ± 6.5 ka), which was previously discussed with regards to mapping, has a paleomagnetic secular variation with implications for the geochronology of another unit mapped near Al-Madinah. Unit bhu has unique paleomagnetic secular variation values that resemble nearby bau with an inclination of 47.0° and declination of 23.2° (Figs. 9A and 11A). Therefore, we infer both bhu and bau erupted within error of each other at 280.4 ± 6.5 ka.

In addition, multiple eruptions are grouped spatially and temporally. One of the most voluminous (∼1.4 km3) and well-exposed (∼138 km2) eruptive groups consists of eight units that are bracketed between bsh at 188.3 ± 4.2 ka and bnu at 129.0 ± 8.2 ka (Fig. 11B). Half of these eight units have age determinations, which help to constrain the remaining undated units. Unit bsh (188.3 ± 4.2 ka) provides a maximum age constraint for the overlying unit bjs. Neither of these units is in direct contact with the other six units in this time interval. Of the remaining six units, undated bru is stratigraphically lowest, with a maximum age constraint provided by the underlying unit of bbd (223.4 ± 8.6 ka). Unit bru is overlain by the undated unit of bra, which is overlain by ban1 at 136.7 ± 5.3 ka and undated bmz, which are both overlain by brm at 134.4 ± 7.5 ka. Unit bnu (129.0 ± 8.2 ka) directly overlies bru, bra, and bmz. Three of the eight units fall within an excursional timeframe. These are bsh (188.3 ± 4.2 ka) and undated bru and bjs. Unit bsh falls into a time interval often referred to as the Iceland Basin excursion (Oda, 2005), while undated bru and bjs have similar paleomagnetic directions (inclinations of 46.3° and 48.5°, declinations of 332.4° and 331.8°, respectively), indicating that they were also erupted during an excursional time interval (Fig. 11B). The excursional nature of these latter two units narrows the time window for when the bru and bjs eruptions occurred to: (1) during the Iceland Basin excursion but after eruption of bsh (188.3 ± 4.2 ka) when the geomagnetic field had sufficient time to move to a different excursional direction; or (2) a subsequent excursion younger than the Iceland Basin event but older than ca. 137 ka. Of the remaining five units, ban1 (136.7 ± 5.3 ka) does not match the remanent magnetization of any other unit, which indicates that it does not cluster temporally with the rest of the group (Fig. 11B). However, undated units bra and bmz, brm at 134.4 ± 7.5 ka, and bnu at 129.0 ± 8.2 ka form a closely spaced array using their average paleomagnetic remanent directions (Fig. 11B). It is therefore suggested that these basaltic units erupted over a relatively narrow time interval of ∼200–500 years based on usual rates of geomagnetic secular variation (Champion and Shoemaker, 1977).

Age determinations indicate that volcanism near the present-day site of Al-Madinah occurred from 1014 ± 14 ka to the historical eruption at 1256 A.D. Only a few mapped units are older than 500 ka, which include bre (1014 ± 14 ka), a reverse polarity scoria cone with no exposed lava flow, bpa (537.1 ± 9.4 ka), and potentially the undated unit bha (Figs. 4 and 9A). Whether the scarcity of exposed units older than 500 ka is due to their burial and concealment by younger eruptives, or if volcanism within our map area has been more active since ca. 500 ka, is currently unknown (Fig. 9B). Ten scoria cones within Figure 4 lack known correlative lava flows, and a similar number of lava flows lack identified scoria cones. Additionally, the thickness of volcanic strata throughout our small map area has yet to be determined through geophysical methods or drilling. Aeromagnetic data collected at greater than kilometer-scale resolution (Blank and Sadek, 1983) was used to create an isopach map of Harrat Rahat in Camp and Roobol (1989), which suggests that most of the area around Al-Madinah is covered by <100 m thickness of mafic lavas, although parts may attain thicknesses up to ∼200 m. While the thickness of lavas at the far-northern margin is regarded as relatively thin (presumed to be a single flow thick) compared to farther south (potentially >300 m thick in areas), new fine-scale geophysical surveys or drill holes are needed to create a more accurate isopach map. Since most of the mafic lava flows are considered to be relatively thin (∼10 m thick typically), there could be a significant thickness of older lava flows buried beneath parts of our map area. This is particularly relevant along the southern and eastern margins of the area depicted in Figure 4, where most units erupted from the main vent axis (Fig. 2). The number of eruptions from northernmost Harrat Rahat could be 30%, or more, greater than exposed and designated units based on orphaned scoria cones and lava flows and the potential thickness of buried volcanic strata. Such potentially older concealed eruptive products would likely have vented >500 ka and consist mainly of alkalic or tholeiitic basalts based on a general lack of more evolved hawaiite or mugearite compositions throughout our study area before that time (Fig. 9B). Runge et al. (2014) suggest that only one-third of vents currently remain exposed within Harrat Rahat, but this number seems unlikely within our smaller map area because it is the young leading edge of volcanic field where much of the volcanic strata is a single lava flow deep, and it is predominantly away from the main volcanic axis.

Integrated paleomagnetic data and geochronology demonstrate that some of these eruptions were temporally clustered. Ages of units bhg (376.1 ± 5.7 ka) and bpr (372.5 ± 9.8 ka) coincide within measured uncertainty, and their mean remanent magnetic directions indicate that they erupted no more than a few centuries apart. Although these units crop out only ∼14 km apart, their aerial distributions along the east (bpr) and west (bhg) flanks of our map area make it improbable that they vented from the same location. Furthermore, differences in major-oxide and trace-element compositions and petrography indicate that these were separate, distinct magma batches, which erupted within a short time interval and relatively close proximity. Additionally, some eruptions are clustered in both space and time. Units bhu and bau are both interpreted to have erupted over a short duration at 280.4 ± 6.5 ka, based on nearly identical and unique paleomagnetic variations and chemical similarities. We interpret these attributes to represent derivation from the same mixed magma batch, of which a small volume that formed bau is inferred to have split off from the larger volume bhu magmatic dike during ascent, and erupted ∼2 km away.

Mafic volcanic fields are generally constructed from multiple adjacent or overlapping monogenetic eruptions, wherein a single magma batch erupted once from each vent (Wood, 1979). While each vent within our map area is interpreted to have erupted only once, some encompass a variety of petrographic and geochemical characteristics that are indicative of multiple magma sources having mixed prior to, or during, eruption. Other than bhu and bau, units such as bla and bdw also consist of lava flows and scoria cones with diverse petrographic and geochemical traits. In particular, bla and bdw are more chemically heterogeneous than the other mafic eruptives. Each eruptive unit for which magma mixing is inferred has two clear end-member compositions, and typically a third hybrid composition between these end members (Fig. 8; see Camp et al. [1987] and Murcia et al. [2017] for similar geochemical plots of bla). End members are typically alkalic and tholeiitic basalts. This is most thoroughly documented for the 1256 A.D. bla eruption as presented in Camp et al. (1987), for which 59 samples were analyzed via XRF along the entire 23 km length of the lava flow. Besides significant geochemical diversity, compositions from both Camp et al. (1987) and our study yield linear trends for incompatible-element ratios that indicate mixing, and contain heterogeneous phenocryst assemblages of olivine ± plagioclase with disequilibrium textures.

Other researchers reject magma mixing, and instead attribute geochemical and petrographic differences to the tapping of a vertically differentiating magma column (Murcia et al., 2017). However, Moufti et al. (2012) showed that trace-element plots (i.e., Nd/Sm versus Nd and La/Yb versus La) have convex trends that argue for magmas from different source regions or that have undergone variable degrees of fractionation and mixing. If crystal fractionation was the exclusive process occurring, such trace-element plots would display linear trends (Hofmann and Feigenson, 1983; Weinstein et al., 2006). Magma mixing has been implicated at small-volume mafic volcanic fields within western Saudi Arabia and around the world using geochemical and petrographic data (e.g., Camp et al., 1987, 1991, 1992; Camp and Roobol, 1989; Shaw et al., 2003; Clynne and Muffler, 2010; McGee et al., 2012, 2013; Moufti et al., 2012). More research is required to better understand the end-member components and behavior of magma mixing, such as why some units display magma mixing, whereas other temporally and geographically clustered units appear to be individual, distinct magma batches.

Vent Alignments and Magma Ascent

Volcanic eruptions are the manifestation of dike ascent through the crust, and while eruptions may occur as single or temporally clustered events, it is unknown how tectonism influences diking events. Such questions are important for understanding where and when future eruptions may occur. Eruptive vents are aligned in many volcanic fields, which is typically inferred as being related to coupled magmatic, volcanic, and tectonic processes (e.g., Connor and Conway, 2000; Valentine and Krogh, 2006; Gaffney et al., 2007; Muffler et al., 2011; Runge et al., 2016). Tectonic extension has been invoked for Harrat Rahat (Camp and Roobol, 1989, 1992), although only a single fault (trending N10°W) with ≤10 m displacement has been mapped within the northern third of Harrat Rahat, and inspection of high-resolution lidar shows no offsets in Quaternary alluvial deposits (unit Qal in Fig. 4) adjacent to the volcanic field. The near absence of young faults contrasts considerably with Camp and Roobol (1989), who depicted faults and/or fractures cutting Harrat Rahat in a variety of orientations. However, upon closer inspection both in the field and using lidar, aerial, and satellite imagery, almost all of those proposed faults or fractures coincide with drainage channels enhanced by erosion along lava-flow contacts.

Most notable within the central-eastern part of northern Harrat Rahat is the high-standing alignment of vents (scoria cones, domes, and craters) of the main vent axis for the northern part of Harrat Rahat (Fig. 2). Scoria cones within our map area, and throughout Harrat Rahat, are broadly aligned with both the main vent axis and Red Sea rift with orientations of about N30°W (El Difrawy et al., 2013). Figure 4 encompasses both elongate scoria cones and aligned craters and cones from individual eruptions, which are the result of common fissure eruptions that can have surface expressions >2 km long (i.e., bla). These vents are oriented between N16°W to N33°W, with an average alignment of N20°W, including the seven cones and craters of the 1256 A.D. bla eruption; these cones and craters have an orientation of N30°W. These vent alignments are used to infer the presence and orientation of subsurface structures that dikes utilized during ascent through the crust. Due to the absence of young faulting within Harrat Rahat, these aligned vents are proposed as an indicator of the local crustal stress regime. This is supported by vent alignments for individual eruptions, as well as the main vent axis, being broadly parallel to faulting within surrounding Tertiary basalts and Precambrian metasediments and metavolcanics. Camp and Roobol (1992) invoke left-lateral movement across the Makkah-Madinah-Nafud volcanic line to explain the right-stepping en echelon pattern of main vent axes throughout Harrat Rahat. Such an explanation could be invoked on a smaller scale for individual vents themselves throughout Harrat Rahat, Khaybar, Ithnayn, and Kura. Future geodetic surveys could confirm this left-lateral movement; however, counterclockwise rotation of the Arabian plate as a result of rifting along the Red Sea and Gulf of Aden would allow for such movement across these volcanic fields and an en echelon pattern of faults and fractures to exist.

Transporting small-volume magmas through the crust via the development and propagation of fractures is well known, but the state of the crust is important in influencing propagation and dike ascent (e.g., Shaw, 1980; Lister, 1990; Lister and Kerr, 1991). Rapid transport would be required to prevent stalling and solidification within the thick (∼40 km) Precambrian crust beneath Harrat Rahat. Buoyancy-driven magmatic fracturing and a stress field with a positive gradient would assist magma ascent to the surface without stalling (Takada, 1989; Spence and Turcotte, 1990). Measurements of volatile concentrations in western Arabian plate eruptives could inform our assessment of buoyancy-driven ascent, but Hawaiian and Strombolian eruptions of mafic units indicate that volatiles are present. For example, the 1256 A.D. bla eruption underwent fire fountaining to heights of ∼500–1000 m (Kawabata et al., 2015). In some cases, groundwater played a major role in eruption explosivity, because some units include phreatomagmatic eruption products (e.g., Murcia et al., 2015). Despite volatile concentrations being unknown, increases are likely related to fractionation, which may relate to a number of factors (i.e., degree of fractionation, ascent rate, etc.). Volatile driven explosive eruptions are affirmed by clast vesicularity, size, and density for mafic vents throughout our map area. While beyond the scope of this paper, episodes of volcanic activity are likely tied to stress changes in the crust, such as periods of elevated extension or a lithostatic response due to Afar plume-related erosion at the base of the lithosphere (e.g., Elkins-Tanton, 2007; Wang et al., 2015; Konrad et al., 2016). An alternative model could involve forceful dike injection resulting from high magma overpressures (e.g., Rubin, 1995; Traversa et al., 2010), which would account for the lack of regional northward-trending Cenozoic crustal stress indicators beyond the limits of the volcanic fields associated within the Makkah-Madinah-Nafud volcanic line.

Estimating magma ascent rates also has implications for forecasting future hazards. The suspension of xenoliths in erupted magmas and chemical diffusion profiles along xenolith margins have been used in other volcanic settings (e.g., Sparks et al., 2006; Rutherford, 2008); however, xenoliths are essentially nonexistent in Harrat Rahat eruptives. Despite a lack of mantle or crustal xenoliths, based on the aforementioned models, Murcia et al. (2017) suggested that magmas could ascend from ∼100 km depth in a time interval as short as ∼12 days to 5 h. Yet, a 1999 seismic swarm interpreted as magma movement within the lower to middle crust below Harrat Rahat occurred over 17 days (Abdelwahed et al., 2016), while a seismic swarm at Harrat Lunayyir (Fig. 1B) was modeled as a dike intrusion propagating to within ∼2 km of the surface over a ∼3 month period in 2009 (Pallister et al., 2010). Geochemical data indicate that Harrat Rahat mafic eruptives are not primary basalts (i.e., low Mg#, Ni, olivine Fo contents), and we therefore infer that these basalts have spent some residence time in or at the base of the crust. Crystallization differentiation in the immediately subcrustal lithospheric mantle would account for minimal assimilation of Precambrian crust (based on isotopic data in Moufti et al., 2012) at conditions appropriate for clinopyroxene phenocrysts to have crystallized and separated from basaltic magmas.

The northernmost part of the Harrat Rahat volcanic field consists of small-volume, mafic lava flows and their source vent scoria cone and crater complexes with compositions ranging from alkalic and tholeiitic basalt to hawaiite and mugearite. In total, 32 Quaternary mafic eruptive units have been mapped underlying and surrounding the city of Al-Madinah based on integrated field relationships, geochemistry, paleomagnetism, age determinations, digital elevation data, and aerial and satellite imagery. In particular, both geochemical and paleomagnetic analyses confirm the validity of the contacts for the 32 mapped units. Although mapped lava flows range from <1–23-km-long, typical lava flows are roughly 10–15 km long (average of 13.5 km), ∼1–3 km wide, and ∼10 m thick. Geochemical variations indicate fractionation, predominantly of clinopyroxene, olivine (± chromian spinel inclusions), and plagioclase. Our 40Ar/39Ar radiometric and 36Cl cosmogenic surface-exposure ages fall within the documented stratigraphic superposition, with eruptives ranging in age from at least ca. 1 Ma to the most recent Harrat Rahat eruption in 1256 A.D. (unit bla). Contrary to prior investigations, the only confirmed Holocene eruption near Al-Madinah produced the vent complex and lava flow of unit bla, with the next youngest unit being the four small scoria cones in the southwestern suburbs that erupted ∼13 ka (unit bdu). Major times of lava-flow emplacement close to, or within, the present-day site of Al-Madinah were during ca. 400–340 ka and ca. 180–100 ka. Radiometric age determinations, coupled with paleomagnetic analyses, indicate that several eruptions were clustered in relatively short timeframes. In particular, bhg at 376.1 ± 5.7 ka and bpr at 372.5 ± 9.8 ka are interpreted to have erupted only a few hundred years apart, while bhu and bau are both interpreted to have erupted at 280.4 ± 6.5 ka from the same magma batch that split during ascent. Similarly, four mafic eruptives (units bra, bmz, brm, and bnu) are proposed to have erupted over at most a few centuries based on stratigraphic superposition, overlapping age determinations, and similarities in paleomagnetic analyses that can be explained through normal rates of geomagnetic secular variation. The relatively primitive nature of these mafic eruptives argues against appreciable assimilation of crustal rocks, partial melts, or the existence of a sizeable, long-lived crustal magmatic reservoir, although the mixing of magmas prior to or during eruption occasionally occurs. These magmas are inferred to follow middle- to upper-crustal structures during ascent to the surface, as evidenced by northwest- to north-trending alignments of scoria cones, craters, and fissure vents. This is consistent with an extensional system along the Makkah-Madinah-Nafud volcanic line, which overlies asthenospheric flow of the Afar plume, although geomorphic evidence is lacking for ongoing upper-crustal extension. As such, forced dike injection due to high magma overpressures represents an alternative model for magma transport to the surface.

This research was funded by the Saudi Geological Survey through a technical cooperative agreement between the Saudi Geological Survey and the U.S. Geological Survey. Thomas Sisson, Andrew Calvert, Joel Robinson, Tim Orr, Fawaz Muquyyim, Mahmod Ashur, and many others are thanked for help with field work. The hard work of James Saburomaru, Dean Miller, Katie Sullivan, and Brandon Swanson in getting samples prepared for 40Ar/39Ar dating experiments is greatly appreciated. We acknowledge Andrew Calvert and Thomas Sisson for internal U.S. Geological Survey reviews. Thoughtful reviews by two anonymous reviewers, associate editor Jan Lindsay, and editorial handling by Shan de Silva greatly improved this paper. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the United States government.

1Supplemental File. Detailed X-ray fluorescence (XRF) and inductively coupled plasma–mass spectrometry (ICP-MS) methods and data (Table S1), detailed 40Ar/39Ar radiometric and 36Cl cosmogenic surface-exposure dating methods, complete 40Ar/39Ar age profiles (Fig. S1). Please visit http://doi.org/10.1130/GES01625.S1 or the full-text article on www.gsapubs.org to view the Supplemental File.
Science Editor: Shanaka de Silva.
Associate Editor: Jan Marie Lindsay.
Gold Open Access: This paper is published under the terms of the CC-BY-NC license.

Supplementary data