Mantle earthquakes in recently thinned Neoproterozoic lithosphere: Harrat Lunayyir, Saudi Arabia

We use earthquake geothermometry, measured heat flow, and structural constraints from P-wave receiver functions to model the thermal evolution of the lithosphere beneath Harrat Lunayyir. We suggest that the lithosphere thinned to its present 60-km thickness in a second stage of lithospheric thinning at 15–12 Ma following initial Red Sea extension at ~27 Ma. Harrat Lunayyir is an active volcanic field located in the Arabian Shield > 150 km east of the Red Sea rift axis. In the lithospheric mantle beneath Harrat Lunayyir we locate 64 high-frequency earthquakes at depths of 42–48 km, all with mL ≤ 2.5. These brittle-failure earthquakes must have nucleated at relatively low temperatures, based upon global maximum nucleation depths and temperaturedependent-deformation experimental results. Therefore, the mantle earthquakes show that the upper-mantle lithosphere is not in thermal equilibrium with the shallow (60 km) underlying asthenosphere. Our thermal modeling indicates that the lithosphere beneath Harrat Lunayyir thinned to its current 60-km thickness at 12 ± 2 Ma, as constrained by thermal modeling of: (1) surface heat-flow; (2) the depth to the mid-crustal brittle-ductile transition, and (3) the depth to the upper-mantle brittle-ductile transition.


Introduction
Earthquakes within the continental mantle lithosphere are sufficiently rare as to excite intensive study. These events are mostly found within zones of continental convergence or subduction (Chen & Molnar, 1983;Sloan & Jackson 2012;Schulte-Pelkum et al., 2019), beneath non-magmatic rift zones (Yang & Chen, 2010), and possibly also plate-boundary strike-slip faults (Inbal et al., 2016). The seismicity of the Arabian plate is highly concentrated along its borders (Youssef, 2015), and the seismic quiescence of the Arabian interior conforms to the very definition of a craton (Mooney et al., 2012). Precambrian shields typically have sufficiently thick and cold lithosphere that mantle earthquakes could occur, but their exceedingly low strainrates, including those within Arabia (Reilinger & McClusky, 2011), cause mantle earthquakes to be particularly rare beneath shields (Frohlich et al., 2015). Volcano-seismicity, with mantle and lower-crustal earthquakes likely caused by early-stage pooling and mixing of asthenospheric melt at the base of or within the lower-crust (McCausland et al., 2017), should also be rare in cratons, including Arabia, given their low heat-flow (Gosnold, 2011). However, here we report 64 high-frequency (brittle-failure) earthquakes within lithospheric mantle beneath the western margin of the Arabian Shield.  (Gosnold, 2011) and are labeled as follows: 1 = SAU002 (34 mW/m 2 ), 2 = SAU003 (43 mW/m 2 ), 3 = SAU004 (44 mW/m 2 ), 4 = SAU006 (50 mW/m 2 ). Colored circles indicate estimates of depth to LAB from Hansen et al. (2007). MMN = Makkah-Madinah-Nafud line.
Understanding of Red Sea rift history has evolved beyond early 'two-stage' models for the Red Sea that include two stages of ocean spreading (at 41-34 Ma and at 5-4 Ma; Girdler & Styles, 1974) or two pulses of 'initial' or 'early stage' rifting (at~34 Ma and at 21-25 Ma; Omar & Steckler, 1995) to widespread acceptance that the first formation of Red Sea rift basins was in Late Oligocene time   (Bosworth, 2015). Our two-stage rifting model distinguishes the Late Oligocene initiation from a mid-Miocene second stage of major lithospheric thinning culminating at 12 ± 2 Ma and has much in common with a recent theoretical model that identifies two-stage evolution as a common feature of successful rifts (Brune et al., 2016). In this latter model, lithosphere that is initially cold and thick, hence strong, begins to rift slowly with full plate-separation rates < 10 mm/yr. As the lithosphere is extended further, it is thinned and heated from below, leading to a nonlinear decrease in lithospheric strength. This thinning and weakening -if it can outpace conductive cooling of the lithosphere -initiates a second, much faster phase of rifting. This second stage of rifting has a profound effect upon the thermal state of rifted margins and hence the depths of their brittle-ductile transitions.
Harrat Lunayyir is a volcanic field along the western edge of the Arabian Shield, east of the Red Sea rift margin (Figs. 1 & 2), beneath which we re-located 64 earthquakes that occurred within the lithospheric mantle during 2014. We assume heat flow within the interior of the Arabian shield at~30-40 mW/m 2 is representative of thermal conditions prior to rifting, whereas a measurement near Harrat Lunayyir is~50 mW/m 2 (Gosnold, 2011). Here we model the thermal evolution of the lithosphere, using constraints provided by depths of earthquake nucleation and measured surface heat flow, to bound the time when the Red Sea rift transitioned from the first slower stage of rifting to the second faster stage.  (Pallister et al., 2010), the local seismometer array, upper-crustal seismicity (depths color-coded red to blue), and mantle seismicity from 2014 (depths color-coded purple to green). The arrow points to the easternmost station of the array (TRAS), 30 km beyond the map. The black box bordering Harrat Lunayyir outlines the extent of our common conversion point cube (Fig. 3). Direction lines within the circle around the mantle earthquakes show the strike of the event cluster (N103°E) and its dip direction (N193°E). The cluster has a dip of 43°from the vertical.

Geologic/Geophysical Framework
The Arabian-Nubian Shield was assembled via repeated terrane accretion during the Neoproterozoic (Stoeser & Camp, 1985). After 500 Myr of relative stability the shield began to break apart due to impingement of the Afar Plume at~30 Ma, which triggered extension first in the Gulf of Aden and then in the Red Sea by 25 Ma (Stern & Johnson, 2010). Initiation of shear on the Dead Sea transform in the middle Miocene and penecontemporaneous initiation of full ocean spreading in the Gulf of Aden facilitated a second phase of Red Sea rifting with rotation of extension from rift-normal (WSW-ENE) to Dead Sea transform-parallel (SSW-NNE) at 11 ± 2 Ma (Bosworth et al., 2015) as constrained geodetically (Reilinger & McClusky, 2011).
Seafloor spreading in the southern Red Sea began at 5 Ma. Spreading has propagated northward to~20°N, but has not yet reached the latitude of Harrat Lunayyir (25°N) (Bosworth et al., 2015). Camp and Roobol (1992) report two distinct phases of volcanism within the Arabian plate. From 30-20 Ma tholeiitic-to-transitional lavas were emplaced along rift-parallel dikes extending the entire length of the Red Sea rift flank, contemporaneous with the impingement of the Afar plume on the base of the African-Arabian plate. From 12 Ma to the present, transitional-to-strongly-alkalic lavas have been emplaced along a N-S trend known as the Makkah-Madinah-Nafud (MMN) line of volcanism forming the younger harrats ( Fig. 1) (Camp & Roobol, 1992;Stern & Johnson, 2010). Two stages of rifting and volcanism are also potentially reflected within the uplift history of the Arabian rift flank, with early uplift of a broad region of the Arabian-Nubian Shield progressively focusing to a narrow region around the Red Sea rift (Szymanski et al., 2016). Mooney et al. (1985) analyzed seismic wide-angle reflection/refraction travel-times along a line through heat-flow measurement points 1, 2, 3 ( Fig. 1) to the Red Sea and constrained the Arabian Shield upper-crustal P-wave velocity of 6.3 km/s and mid-to-lower crustal velocities of 6.5-6.7 km/s, with the Moho~40 km depth. Not far below the Moho, both body-wave (e.g. Koulakov et al., 2016) and Rayleigh-wave (e.g. Yao et al., 2017) tomography show a strong decrease in S-wave velocity along the entire western margin of the Arabian craton, below~60 km beneath Harrat Lunayyir and gradually deepening toward the plate interior. Comparison with Swave receiver function converters (Hansen et al., 2007;Fig. 1) suggests this S-wave velocity decrease marks the lithosphere-asthenosphere boundary (LAB).
Volcanism within Harrat Lunayyir began at~600 ka sourced by melt segregation within the asthenosphere at~60-75 km depth, and since then eruption rates have decreased and depth of melt segregation may have decreased (Duncan & Al-Amri, 2013). Harrat Lunayyir last erupted during the 10 th century C.E. In 2009, a roughly N-S striking dike (Fig. 2) of~0.13 km 3 reached within 1-2 km of the surface, and the associated >30,000 earthquakes, all with mb ≤ 5.7, caused evacuation of the nearby town of Al-Ays (Pallister et al., 2010). This volcano-tectonic crisis led to densification of the existing seismometer array from three sites before 2009 to 18 sites by August 2010, which provide the data for our study.

P-Wave Receiver Functions
The Harrat Lunayyir seismic array now consists of 18 Trillium 120 broadband stations (Table 1) and spans a region~50  50 km (Fig. 2). We used data from this array to create P-wave receiver functions (PRFs) from 857 earthquakes recorded during 2008-2015 with epicentral distances of 30-90°, largely from backazimuths of 020-120°and with moment magnitudes of 5.5-8.3 (Fig. 3). P-wave receiver functions (Langston, 1979) are created for each station using teleseismic events recorded from 2009-2014. The three-component seismograms were first rotated to the Radial-Transverse-Vertical (RTZ) coordinate system to isolate the first arrival on the vertical component and the P-to S-wave conversion on the radial component. Arrival times for every station-event pair were initially picked using a short-term average (STA) over long-term average (LTA) algorithm based on the method of Earle and Shearer (1994). These initial arrival times were then manually refined to correct for picking errors, as well as visually inspected to ensure that data with no clear arrivals were excluded from the analysis.
To calculate each receiver function we took 100 s of seismic data, filtered from 0.2 to 1.5 Hz. The traces begin 20 s before the first P-wave arrival and include 80 s of post-arrival data. The receiver functions were constructed by deconvolving the vertical (Z) signal from the radial (R) trace using the iterative-time domain approach of Ligorria and Ammon (1999) as implemented in the codes of Hermann (2013). Receiver functions were visually evaluated for consistency by comparing each event at a station to all other events at the same station. We then performed common conversion-point stacking of the receiver functions in order to gain a better estimate of the three-dimensional structure of the lithosphere using the FuncLab program (Eager & Fouch, 2012). Each stack was bootstrap-resampled by 150% and a minimum threshold value for plotting each voxel was set equal to the standard deviation of the stack at each voxel. Our common conversion-point "cube" (Fig. 3) shows a clear, coherent and continuous positive polarity discontinuity (i.e., an upward decrease in seismic impedance for the P-to-S conversion) at 38 km below the surface. The average half-width of the pulse at half-maximum amplitude provides our uncertainty estimate of ± 2 km. This positive polarity ( Fig. 3) represents the Moho, in agreement with previous crustal thickness estimates (Hansen et al., 2007;Park et al., 2008;Mechie et al., 2013;Tang et al., 2016). At 60 km we observe a seismic converter with opposite polarity to that of the Moho (Fig. 3) that we interpret as the lithosphereasthenosphere boundary (LAB), also consistent with previous lower-resolution estimates (Hansen et al., 2007;Park et al., 2008;Koulakov et al., 2016). The average half-width at half-maximum of the LAB signal is also~2 km, but due to its less-focused appearance we use a more conservative uncertainty estimate of ± 5 km. It is uncommon to image the LAB using P-wave receiver functions due to the interference of crustal multiples (Kind et al., 2012), but here we benefit from dense station coverage, numerous teleseismic events, and a mantle-lid that is significantly thinner than the crust. In stark contrast to the 60-km thickness beneath Harrat Lunayyir, a lithospheric thickness of 100-150 km is reported at the transition from the shield to the platform beneath central Saudi Arabia (Hansen et al., 2007) and even thicker lithosphere (≥ 200 km) beneath the eastern Arabian Peninsula (Koulakov et al., 2016).  Figure 2. The color-scale represents the amplitude of stacked conversions (P-to S-waves) with respect to the incident P-waves. Negative amplitudes (blue) indicate a wave-speed decrease with depth and positive amplitudes (red) indicate a wave-speed increase with depth. b) The azimuthal distribution of teleseismic events used to generate Pwave receiver functions. The gold star is centered on the study area.

Earthquake Spatio-temporal Patterns
The Saudi Geological Survey (SGS) provided an event-catalog for 2009-2014 ( Fig. 4) and continuous digital seismic data for the entire calendar year 2014. The 2014 event-catalog contains 72 earthquakes with unusual source depths of 40-50 km beneath Harrat Lunayyir and with local magnitudes of 0.2-2.5 (Fig. 5). We re-located the mantle earthquakes using hypoinverse (Klein, 2002) to improve absolute locations and HypoDD (Waldhauser, 2001) to reduce relative errors between events, using the IASP91 velocity model (Table 2; Kennett, 1991). Our relocations confirm that 64 of these catalog events are indeed below the 38-km deep Moho, all at depths of 42-48 km and all located~30 km southeast of the 2009 dike ( Fig. 2 & 5; Table 3). Reported relative depth uncertainties from HypoDD (Waldhauser, 2001), using the SVD method, are ± 100 m. Absolute uncertainties depend on array geometry. Our preferred earthquake depths minimize the observed S-minus Pwave travel-time residuals with a median RMS error of < 0.23 s ( Fig. 6 & 7) across all 64 sub-Moho events. Because there exists only one station east of these earthquakes (TRAS, Fig. 6 & 7), we calculate the RMS of Sminus-P residuals (i.e., observed-minus-calculated travel-times) from all stations for each event, on a 1-km grid across a 1000 km 2 N-S swath encompassing the hypocenters and at depths every 0.33 km from the surface tõ 70 km (i.e. 210,000 trial locations) (Fig. 6) and across a W-E swath (Fig. 7). A histogram of the number of source locations yielding a given RMS S-minus-P residual, plotted against the difference in depth between the tested source location and our preferred depth ( Fig. 6 & 7), shows a narrow minimum in the distribution of residuals, indicating that our better constrained earthquakes have depth uncertainties < 5 km. Hence, we confirm that all 64 events are likely in the mantle.     We employ one final test to demonstrate that the earthquakes are indeed within the mantle lithosphere. The ability to accurately resolve the hypocentral location of an earthquake is strongly affected by the chosen velocity model, with incorrect models providing spurious results. Schulte-Pelkum et al. (2019) developed an elegant solution to this problem: comparing travel-time differences between S-and P-wave arrivals from local earthquakes with the equivalent travel-time difference from P-wave receiver functions at a nearby station, with both travel-time differences converted to vertical incidence. We used the station closest to the mantle earthquakes (LNYS, Fig. 2 and Table 1) as it is only~7.5 km SW of the centroid of the mantle earthquakes. The comparison between PRFs and earthquake travel time differences in shown in figure 8. The energy from Pto-S converted waves at the Moho interface arrives at~4 seconds, from the LAB at~5.9 seconds, and the S minus P travel times for the mantle earthquakes span 5-5.5 seconds, clearly between the two interfaces. As noted by Schulte-Pelkum et al. (2019), the comparison between S minus P travel-time differences is rather insensitive to changes in velocity model so this technique clearly verifies our events are depper than the Moho. In figure 8 the mantle earthquakes are closer in S-minus-P time to the LAB than to the Moho, but this is an artifact due to the large change in wavespeed across the Moho. Figure 8 shows the relative, not the absolute depth of the mantle earthquakes.    Fig. 2) and receiver function piercing points (transparent white diamonds) calculated at a depth of 40 km through the IASP91 reference velocity model. b) The azimuthal distribution of teleseismic events used to generate P-wave receiver functions, isolated to the same back-azimuth range as the mantle earthquakes to seismic station LNYS (117 in total). The gold star is centered on the study area. c) Comparison of vertical S minus P travel times from P-wave receiver functions (background) to mantle earthquakes (transparent gray circles, 59 in total), corrected to vertical incidence using a Vp = 6.5 km/s and assuming a Poisson solid.
Also during 2014,~6,000 earthquakes occurred at depths between 0-20 km, with the vast majority of events spatially associated with the 2009 dike ( Fig. 2 & 5). There are three earthquakes in the 2014 SGS catalog with lower-crustal depths of 20-38 km. The waveforms for all three "lower-crustal" events have emergent/noisy first arrivals, making precise locations difficult to determine. We conclude that the lower-crust at > 20 km depth beneath Harrat Lunayyir was essentially aseismic in 2014. Figure 9. a) Three-component seismogram recorded by station LNY01, instrument response removed, of a mL = 1.1 earthquake that has a source-depth of 45.4 km, high-pass filtered above 1 Hz. b) Spectrogram of the vertical (Z) seismogram showing peak P-wave energy at~30 Hz and S-wave energy between 10-20 Hz. Nyquist frequency is 50 Hz.
The mantle earthquakes beneath Harrat Lunayyir are dominated by high-frequency (> 5 Hz) energy, with a peak at~30 Hz (Fig. 9). They have clear, impulsive P-and S-wave arrivals and match the characteristics of "brittle-failure" earthquakes (Hill & Prejean, 2005). Brittle-failure earthquakes close to sites of active magmatism are commonly referred to as volcano-tectonic because they are indistinguishable from tectonic earthquakes despite occurring around volcanoes. Brittle-failure earthquakes contrast both with long-period (LP) volcano seismicity (0.5-5 Hz), caused by resonance within fluid-filled cracks, and with very long-period (VLP) seismicity (0.01-0.5 Hz), associated with the unsteady transport of magma and gases through conduits (Chouet & Matoza, 2013). High-frequency, LP and VLP crustal earthquakes all occurred during the 2009 diking episode (Pallister et al., 2010), but none of the 2014 mantle earthquakes we examined have LP or VLP characteristics. The 2014 Harrat Lunayyir earthquakes within the mantle lithosphere follow a three-part temporal evolution (Fig. 10b & d) along a plane striking~N103°E and dipping~43°S (Fig. 5). First, between January and early May 2014, events with magnitudes between 0.5-1.5 occurred at a mean depth of~45 km. Then, during an eight-day window in May 2014,~50% of the total known earthquakes occurred, 0.6 ≤ mL ≤ 1.5, and quickly migrated to the shallowest observed depth of~42 km. The subsequent month-long quiescence was broken in June 2014 when the deepest sequence began at~46 km, increasing to 48 km, and encompassing the entire magnitude range (0.2-2.5) of the 2014 mantle earthquakes.
This three-part temporal trend does not correlate with the upper-crustal seismicity (Fig. 9d), which is roughly constant throughout the entire year. The temporal evolution of the earthquake locations matches "stage 1" volcano seismicity caused by magma-driven fracturing and pooling between the LAB and Moho (McCausland et al., 2017). Similar patterns of seismicity are present in the SGS catalog between August 2010 and the end of 2013 (Fig. 10a & c), indicating that the 2014 sequence is not an isolated occurrence and that there are likely numerous pulses of magma intruding close to the base of the crust.

Derivation of Finite Differences Thermal Modeling Scheme
We used a one-dimensional finite-differences code to model the thermal evolution of the lithosphere as it is eroded from below, following Turcotte and Schubert (2002) and Recktenwald (2004). We start with the one-dimensional time-dependent heat-conduction equation (Turcotte & Schubert, 2002, eqn. 4.68), with the inclusion of radiogenic heat production: where = 㘠 is the thermal diffusivity and 㔵 is the volumetric heat production of the medium (for parameter definitions see Table 4). If we assume that the concentration of heat-producing elements decays exponentially as a function of depth from the surface with a characteristic decay constant hr, we can write the volumetric heat production as: (2) into equation (1) we obtain: where = 0 is a constant related to the medium. This is the one-dimensional heat-flow equation with internal heat production that we next solve using the centered-space finite-differences approximation for the term involving the second spatial derivative of the temperature: where m indicates the spatial node on a discrete mesh. Next, we use the forward-time finite-differences approximation for the left-hand side of equation (3): where n indicates the temporal node on a discrete mesh. In equations (4) and (5), z and t are the distances between spatial and temporal nodes, respectively. Substituting equations (4) and (5) into equation (3) produces: where superscripts denote temporal nodes and subscripts are spatial nodes. Rearranging the terms to solve for the temperature at spatial node m in the next time step, n+1, we get: We simplify the equation by setting: = ∆㠴 ∆㔵 2 ; = 1 − 2 ; = ∆㠴 to obtain: # 8 This is an explicit equation for the temperature at spatial node m and time step n+1 only in terms of the previous time step. The constant r is the stability parameter for the finite differences scheme. The solution is stable on the condition that ≤ 1 2 (Recktenwald, 2004). To obtain an input temperature profile, we use the solution for the steady-state heat-flow equation with internal heat generation (Turcotte & Schubert, 2002, eqn. 4.30): To solve for qm, the reduced heat flow, we apply the boundary condition at the surface (z=0) as shown in equations (10) and (11): − 㔵 = 0 = 0 = + 0 ℎ # 10 = 0 − 0 ℎ # 11 In order to use equation (8) to solve equation (3), we need two boundary conditions: temperatures at the surface (T0) and at the lithosphere-asthenosphere boundary (TLAB at depth zLAB). The surface boundary condition is simply the average temperature at the surface (25°C, Table 4). Mantle temperature depends on the mantle potential temperature, Tp, and the adiabatic temperature gradient within the mantle:

Global Heat Flow Database
We assume that the seismic lithosphere-asthenosphere boundary is coincident with the thermal definition of the lithosphere-asthenosphere boundary as the intersection between the conductive lithospheric thermal profile and the adiabatic mantle-asthenosphere thermal profile. If, as we show in this paper, the lithosphere is not in thermal equilibrium with the underlying asthenosphere due to geologically recent thermal erosion, then we would indeed expect a narrow seismic transition zone (seismic LAB) at the thermal lithosphere-asthenosphere boundary (Aulbach et al., 2017).

Thermal Constraints
The maximum depth at which intra-continental earthquakes nucleate is controlled by strain-rate, lithology and temperature (e.g. Chen & Molnar, 1983). Reilinger and McClusky (2011) report that strain-rates within the Arabian Shield are below their detection threshold (< 10 -15 s -1 ). Due to our interest in the mantle lithosphere, we restrict our analysis to olivine-dominated lithologies. The clear cut-off in earthquakes at 48 ± 2 km depth below Harrat Lunayyir shows that at this depth the modern temperature is at the maximum at which brittle-failure earthquakes can occur, which must be ≤ 900-1000°C as measured from xenoliths at Harrat Kishb (McGuire & Bohannon, 1989). Chen & Molnar (1983) suggest that crustal and mantle seismicity are limited to temperatures less thañ 350 ± 100 and~700 ± 100°C, respectively, thought to represent brittle-ductile transitions in the crust and upper mantle. More recently, a 600°C cut-off in seismicity in the mantle has been estimated for oceans (McKenzie et al., 2005) and continents (Sloan & Jackson, 2012). However, in volcanic areas, due to elevated strain-rates, lower-crustal brittle-failure earthquakes are well known from regions thought to be hotter, 700-750°C (Shelly & Hill 2011;White et al., 2011). Ágústsson & Flóvenz (2005) used geothermal data to constrain the maximum nucleation depths of earthquakes beneath Iceland to temperatures of 600-900°C. Although lithological considerations alone suggest mantle-lithosphere brittle-failure earthquakes should therefore be common in continental volcanic areas, they are only rarely reported. For example, there are no brittle-failure events beneath the Cascade volcanic arc where LP earthquakes are recognized around Moho depth at temperatures in excess of~800°C (Hansen et al., 2016). Beneath Hawaii, tectonic earthquakes thought to be caused by volcano loading and lithospheric flexure are known to occur as deep as 50 km along a long-lived tectonic fault zone (Wolfe et al., 2003). Recently, Inbal et al. (2016) located small earthquakes in southern California along the Newport-Inglewood fault below the Moho where temperatures are~800°C, correlating with geochemical evidence of a fluid pathway from the mantle (Boles et al., 2015). Prieto et al. (2017) suggest that a mantle earthquake beneath the Wyoming craton at 75-km depth occurred at temperatures of 750-850°C due to strain localization in the nominally ductile regime,~100 km away from the nearest volcanic center.
Laboratory experiments to determine rate-and-state friction laws for olivine are consistent with these earthquake observations, suggesting a temperature cut-off of~600°C at strain rates of 10 -15 -10 -14 s -1 for typical earthquakes, with the possibility that low-asperity stresses could extend this temperature to~800°C (Boettcher et al., 2007;King and Marone, 2012). Ohuchi et al. (2017) detected faulting associated with accelerated strain localization and localized heating in very fine-grained shear zones in dry dunite and extrapolate their laboratory observations to imply the possibility of faulting in the mantle lithosphere at temperatures up to 900°C at strain rates as low as 10 -14 s -1 . Thus, earthquake geothermometry suggests temperatures of 800°C, or even 900°C, as the cause of our observed 48-km seismicity cut-off.
There are no xenoliths reported from Harrat Lunayyir, but at Pleistocene Harrat Kishb, 150 km east of the MMN line (Fig. 1), spinel peridotites suggest temperatures of 910-970°C at 13-20 kb (McGuire & Bohannon, 1989), corresponding to the entire depth range from the Moho to the LAB as imaged in this paper. Seismic wave-speeds, generally only a weak indication of temperature, suggest the mantle lithosphere beneath Harrat Lunayyir is at 800 ± 200°C (Tang et al., 2016). Finally, our observation of seismograms with spectrogram peaks > 15 Hz (corresponding to short-wavelength S-waves, < 300 m) (Fig. 6) indicates moderate to high Q in the source region. The sharp increase in seismic attenuation within olivine-dominated rocks at 1000°C (Artemieva et al., 2004) allows us to confidently exclude temperatures ≥ 1000°C at depths shallower than the seismicity cut-off beneath Harrat Lunayyir. We conclude that the ambient temperature at 48 ± 2 km depth beneath Harrat Lunayyir is 850 ± 50°C.
In order to deduce the timing of lithospheric thinning beneath our study area, we assume the thermal structure beneath Harrat Lunayyir prior to any Red Sea rifting or plume impingement matches that in the center of the Arabian craton, which we derive from measured surface heat flow and heat production (measurement #1 ,  Fig 1; Gosnold, 2011) and assuming a mantle potential temperature of 1354°C (as petrologically estimated on the flanks of the MMN line; Camp & Roobol, 1992). We solve for the depth to the LAB numerically and find the original depth to be 113-156 km, in accord with the modern 100-150 km values from S-wave receiver functions at the transition from the shield to the platform (Hansen et al., 2007). We use 1D finite-differences to model the temperature evolution beneath Harrat Lunayyir as and after the lithosphere is thinned to our observed 60 km thickness (Fig. 11). We assume the lithosphere is replaced from below by asthenosphere with an adiabatic gradient of 0.3°C/km, and we include depth-dependent radiogenic heat production and thermal conduction throughout model time and space. Our simple thermal model is insensitive to the duration of thinning so long as the duration is less than~70 Myr, because conduction has insufficient time to impact isotherms shallower than those being eroded from beneath.
Using the 350 ± 50°C and 850 ± 50°C isotherms as cut-off temperatures for seismicity within the uppercrust and mantle lithosphere, respectively, along with a present-day surface heat flow of 50 ± 4 mW/m 2 (measurement #4 , Fig 1; Gosnold, 2011), we can place upper (14 Ma) and lower (10 Ma) bounds on how long ago the LAB must have reached its present, shallow depth of 60 km. Figure 11. Thermal evolution model for the lithosphere beneath Harrat Lunayyir. Negative time is associated with the thinning process and positive time indicates the amount of time that has elapsed since thinning completed. a) Surface heat-flow (SHF) predicted by our model through time. Red line marks 50 mW/m 2 (modern heat-flow) and gray lines mark the ± 4 mW/m 2 bounds (Gosnold, 2011). b) Model run with an initial geothermal profile at -10 Myr, matching current measurements from the center of the craton (i.e. Heat Flow 1, Fig. 1). c) Enlarged view of the timeframe in which the 350 ± 50°C isotherm obtains a depth of 20 km (cut-off depth of the shallow seismicity). d) Enlarged view of the timeframe in which the 850 ± 50°C isotherm obtains a depth of 48 ± 2 km (cut-off depth of mantle seismicity). Thinning of the lithosphere stops with the LAB at 60 km at time t = 0 Myr. Initial lithospheric thickness is the intersection of the conductive geotherm and the mantle adiabat. We thin the lithosphere to 80 km by basal erosion by -1 Myr, and then thin the lithosphere very rapidly from 80 to 60 km from -1 to 0 Myr.

Sensitivity of Thermal Modeling Parameters
We tested the sensitivity of our thermal model to various parameter choices and assumptions. In our model of Figure 11 we used the commonly cited crustal conductivity of 2.51 W/m°C (e.g. Gettings, 1982). Whittington et al. (2009) suggest thermal conductivity values as low as~2 W/m°C in hot lower crust; however, using lower values of thermal conductivity raises the temperatures within the lithospheric mantle, so that we cannot match our observations of brittle-failure earthquakes within the lithospheric mantle unless such earthquakes can nucleate above 900°C. Figure 12. Comparison of thermal evolution models for mantle potential temperature of 1354°C and initial surface heat flow of (a) 30 mW/m 2 and (b) 34 mW/m 2 . Top) Calculated surface heat flow from the model; red line marks 50 mW/m 2 (modern heat flow) and gray lines mark the ± 4 mW/m 2 bounds. Pink overlay marks the timeframe when the modeled heat flow values match the measured values. Middle) Thermal model at depths of 15-30 km. Gray line at 20 km depth is the cut-off depth of the shallow seismicity. Blue overlay marks the timeframe when the base of the crustal seismicity (20 km) is within the temperature range of 350 ± 50°C. Bottom) Thermal model at depths of 40-55 km. Thick gray line is the 48 km cut-off of the mantle seismicity, with the ± 2 km bounds above and below it. Yellow overlay marks the timeframe when the base of the mantle seismicity (48 ± 2 km) is within the temperature range of 850 ± 50°C. The vertical gray bar marks the timeframe where all the constraints overlap, denoting the preferred timing of lithospheric thinning to 60 km.
In all our models the predicted initial surface heat-flow (SHF) (e.g. 42 mW/m 2 in Fig. 11a) exceeds the specified initial SHF (30 mW/m 2 for Fig. 11a) because the predicted SHF is calculated using the near-surface thermal conductivity actually measured at point #4 ( Fig. 1) (ksurf = 3.51 W/m°C) since we are trying to match the heat-flow measured at point #4 (50±4 mW/m 2 ) (M. Gettings, unpublished manuscript;Gosnold, 2011). In contrast the thermal calculations (e.g. Fig. 11b) use a bulk lithospheric thermal conductivity of kavg = 2.51 W/m°C. Thus, the model shown in Figure 11 (identical to that in Fig. 12a) is chosen to have an initial SHF of 30 mW/m 2 for a conductivity of 2.51 W/m°C and hence a near-surface thermal gradient of 12°C/m. We use this near-surface thermal gradient (12°C/m) with ksurf = 3.51 W/m°C to calculate the predicted SHF plotted in Figure 11a (42 mW/m 2 ).
We tested the sensitivity of our thermal model to a larger initial SHF of 34 mW/m 2 (the best estimate of Gettings & Showail (1982) for measurement #1 is 34 ± 4 mW/m 2 ) (Fig. 12b), which predicts SHF of 48 mW/m 2 at 0 Myr model time due to the high shallow thermal conductivity. Because this predicted SHF of 48 mW/m 2 is already within the uncertainty of the observations at measurement #4 (Fig. 12b), observed surface heat flow cannot place a lower limit on the age of lithospheric thinning. Nonetheless, the depths to the two brittle-ductile transitions still provide a strong constraint on the time of final thinning, 10-11 Ma for the initial heat flow of 34 mW/m 2 (Fig. 12b), very close to our previous result of 12-13 Ma for initial heat flow of 30 mW/m 2 (Fig. 12a).
The above results are all calculated for a mantle potential temperature of 1354°C. If we instead use a mantle potential temperature of 1200°C, lower than plausible, thinning must complete by~15 Ma for initial heat flow of 30 mW/m 2 (Fig. 13a), or by~12 Ma for initial heat flow of 34 mW/m 2 (Fig. 13b).
Thus, for the likely range of initial surface heat flow values and the mantle potential temperatures we can limit the time-frame of thinning to approximately 12 ± 2 Ma.
Finally, we consider the consequences for our thermal modeling if we have incorrectly identified the prominent negative polarity (blue) converter in our common conversion-point cube (Fig. 3) as the LAB when in fact it could be a 'mid-lithospheric discontinuity' (MLD) (Aulbach et al., 2017). Beneath Harrat Lunayyir we have strong evidence from both body-wave and surface-wave tomography (e.g. Koulakov et al., 2016;Yao & Mooney, 2017) that our converter is at an appropriate depth to be the LAB. However, it is not impossible that an MLD could be present at 60 km depth, if so, presumably above a deeper LAB invisible to receiver functions. Such an MLD would perhaps represent a live or frozen melt percolation front of carbonated peridotite representing solidus temperatures at 60 km depth in the range 1050-1150°C (that would correspond to equilibrium geotherms of~60-70 mW/m 2 ) (Aulbach et al., 2017;Dasgupta et al., 2013). We therefore ran our thermal model for a "mantle potential temperature" of 1050°C (Fig. 14) to show that even in this unlikely case the model results are not vastly dissimilar to those with our preferred parameters: thinning must complete bỹ 18 Ma for initial heat flow of 30 mW/m 2 (Fig. 14a), or by~14 Ma for initial heat flow of 34 mW/m 2 (Fig. 14b).  In summary, our three independent thermal constraints (thinning must be recent enough that a brittle mantle lid still exists; thinning must be old enough that the crustal brittle-ductile transition has begun to shallow; and just right to match the observed surface heat-flow) provide a temporal constraint that is rather insensitive to a wide range of initial thermal conditions. Figure 11d shows that by using the full range of geothermal models from Gosnold (2011), a temperature of 800°C (900°C) at 48 ± 2 km depth is reached 2-5 Myr (4-11 Myr) after the lithosphere has thinned to 60 km depth. Although this time range is broad (2-11 Myr), it is nonetheless instructive. The youngest volcanotectonic event that can possibly have caused local lithospheric thinning, the onset of Harrat Lunayyir volcanism, is dated at 600 ka (Duncan & Al-Amri, 2013). However, if the lithosphere thinned to 60 km as recently as 600 ka, the temperature at 48 km depth would now be only 600 ± 50°C (geotherm at time = 0.6 Myr, Fig. 11b), and we would expect to see deeper earthquakes within the mantle lithosphere. The surface heat flow would also not have had time to increase from the initial 30 mW/m 2 to 50 mW/m 2 . Hence, the onset of Harrat Lunayyir volcanism at 600 ka was not the cause of local thinning of the lithosphere to the modern LAB depth.

Results and Discussion
The oldest magmatic event likely to have thinned the cratonic lithosphere was the impingement of Afar plume-head asthenosphere at~30 Ma and the initial extension of the Red Sea by 25 Ma (Stern & Johnson, 2010;Bosworth, 2015). However, if the lithosphere beneath Harrat Lunayyir had thinned to its present 60 km as early as 25 Ma, the temperature at 48 km depth would now be 1050 ± 50°C (geotherm at time = 25 Myr, Fig. 11b), and we would not have brittle-failure earthquakes at that depth nor most likely anywhere in the mantle lid. Furthermore, we would expect high seismic attenuation, inconsistent with our observation of high-frequency earthquakes, as well as significantly elevated surface heat flow values of~70 mW/m 2 , far higher than observed (Fig. 11a). Hence, early Red Sea extension cannot have thinned the Arabian lithosphere to its present 60-km thickness beneath Harrat Lunayyir.
The implication of the brittle-failure mantle earthquakes is that a second phase of lithospheric thinning, as well as uplift and likely extension, must have taken place in the Middle and/or Late Miocene. The younger harrat volcanism along the MMN line began at~12 Ma (Camp & Roobol, 1992) including other harrats at the latitude of Harrat Lunayyir (e.g., Harrat Kura, 11 Ma; Stern & Johnson, 2010). Apatite and zircon (U-Th)/He ages from the rift margin across Harrats Lunayyir and Khaybar (Fig. 1) suggest renewed extension between 15-12 Ma (Szymanski et al., 2016). If the Arabian lithosphere thinned to 60 km at 12 Ma, then the coldest of our initial geothermal profiles would only now have reached the highest permissible temperature (900°C) at the shallowest permitted seismic cut-off depth (46 km). Our observed earthquake cut-off depth in the mantle is only defined by a relatively short duration of recording; however, if we had observed deeper earthquakes, it would require more recent thinning of the lithosphere. Our choice of a limiting temperature of 850 ± 50°C is higher than the 600-700°C that is often considered the cut-off of brittle behavior in the upper mantle, but using the lower temperatures would only permit a far more recent thinning of the lithosphere. Hence both scenarios (deeper earthquakes; lower temperature limit) would imply that the lithosphere is even more out of equilibrium and would more strongly support two-phase rifting.
The cut-off in crustal seismicity at 20 km depth is far less diagnostic of the timing of the completion of lithospheric thinning because the 350 ± 50°C isotherms in our models are so much flatter than the 850 ± 50°C isotherms (Fig. 11c). Even modest uncertainty in the appropriate cut-off isotherm that arises due to crustal lithologic variability (Magistrale & Zhou, 1996) leads to significant uncertainty in the timing of crustal heating. However, our initial geothermal models, prior to thinning, place the 350 ± 50°C isotherms at~30-45 km depth, clearly deeper than our mid-crustal cut-off in seismicity. Considering the full range of acceptable initial heatflow and heat-production estimates (Gosnold, 2011) permits a crustal cut-off in seismicity at 20 km depth. Hence, we prefer models in which shallowing of the LAB is sufficiently far in the past (before~11 Ma; Fig.  11c) for mid-crustal temperatures to be sufficiently elevated.
Acceptable solutions using the mean initial geothermal parameters for the time at which lithospheric thinning reduced the lithosphere to the observed 60 km thickness must fit our heat flow constraints (8-14 Ma; Fig. 11a), the maximum depth of crustal earthquakes (11-23 Ma; Fig. 11c), and the maximum depth of mantle earthquakes (2-13 Ma; Fig. 11d). Therefore, acceptable solutions must be at least 11 Ma, but can be no older than 13 Ma. Thinning ages < 11 Ma would imply that the crustal seismicity is limited to regions at temperatures < 300°C. (implausible in anhydrous Proterozoic basement; Chen & Molnar, 1983). Thinning ages > 13 Ma would imply that the mantle earthquakes occur at temperatures in excess of 900°C, hotter than previously observed or proposed (Chen & Molnar, 1983;McKenzie et al., 2005). Tests of model sensitivity to various thermal parameters suggest that thinning completed at~12 ± 2 Ma.
Although our modeling cannot formally exclude continuous thinning from 30 Ma (or earlier) to~12 Ma, this 'continuous' scenario offers no ready explanation for the hiatus in magmatism and extension from 20-15 Ma (Camp & Roobol, 1992;Stern & Johnson, 2010). Indeed, a common feature of continental rifts that progress to ocean spreading is that they evolve in two separate phases (Brune et al., 2016), where initially strong lithosphere is slowly thinned and dynamically weakened through rift-induced heating, initiating a second and much faster phase. The second, faster rifting phase usually begins~10 Myr before inception of breakup of the continental lithosphere and persists until complete plate separation is achieved. The second stage is not observed in failed rifts, where conductive cooling strengthens the lithosphere faster than it can be dynamically weakened, leading to rift failure (Brune et al., 2016).
In our preferred scenario the first stage of extension, uplift, and magmatism between 30-20 Ma was contemporaneous with modest lithospheric thinning but failed to rupture the continent and represents the slow phase of rift-evolution. The thermal constraints are, unfortunately, insensitive to the intermediate lithospheric thickness achieved during the first stage given the timescale at which heat diffuses. The change in plate boundary conditions in the middle Miocene (i.e., the onset of full ocean spreading in the Gulf of Aden and of the onset of Dead Sea transform motion) allowed rotation and an increase of the extensional stresses along the Red Sea rift to be more favorable for the second-stage of rifting (Reilinger & McClusky, 2011). Lithospheric thickness beneath Harrat Lunayyir decreased to 60 km between 15-10 Ma, largely by replacement of the lithosphere with asthenosphere, initiating the younger harrat volcanism. This second phase of lithospheric thinning must then have essentially stopped beneath Harrat Lunayyir to preserve the observed gap in seismicity between 48 km and the 60-km deep LAB. We suggest thinning beneath Harrat Lunayyir concluded as rifting became more focused within the Red Sea (cf. Szymanski et al., 2016).

Conclusions
Our seismic imaging locates the Moho beneath Harrat Lunayyir at 38 ± 2 km depth and the LAB at 60 ± 5 km depth, despite being located within a Precambrian shield. Observations of earthquakes within the lithospheric mantle beneath Harrat Lunayyir span the entire time period for which a suitably dense seismic array was active (i.e. since August 2010). During 2014, there were 64 high-frequency "brittle-failure" type earthquakes whose hypocenters, well constrained to the upper 10 km of the lithospheric mantle (depths of 42-48 km), lie on a plane striking~N103°E and dipping~43°S, located roughly 30 km to the southeast of the 2009 dike intrusion. These brittle-failure earthquakes are likely caused by melt migration within the lithospheric mantle. The existence of these brittle-failure earthquakes constrains the temperature at 48 ± 2 km depth to be 850 ± 50°C, which in turn requires the LAB to have reached its present shallow level in middle Miocene time (ca. 12 ± 2 Ma). Our results are consistent with a two-stage model for the evolution of continental rifting, as suggested by Brune et al. (2016) globally and as proposed by Szymanski et al. (2016) for the Red Sea rift. Our identification of high-frequency earthquakes in the shallow lithospheric mantle provides a powerful means of estimating the thermal structure of the lithosphere and consequently the timing of continental rifting in the Red Sea.