{'fraction_non_alphanumeric': 0.04804674033349664, 'fraction_numerical': 0.032678315460962025, 'mean_word_length': 4.409166748454215, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Thickness of the Saudi Arabian Crust

We analyzed P-wave receiver functions from seismic stations covering most of Saudi Arabia to map the thickness of the crust across the Arabian plate. We present an update of crustal-thickness estimates and fill in data gaps for the western shield and the rifted margin at the Red Sea, as well as the eastern Arabian platform. Our application of a conventional H-k stacking algorithm included careful attention to stacking weights, two forms of sedimentary corrections for stations located on the Arabian platform, and additional processing for noisy stations. Average crustal thickness (i.e. depth to Moho below surface) beneath the Red Sea coastal plain (the rift margin) is 29 km, beneath the volcanic harrats is 35 km, the shield (excluding harrats) is 37 km, and the platform is 38 km. Crustal thinning appears not to extend east of the rift escarpment, suggesting uniform extension, that is no broader at depth than at the surface. In contrast to some previous claims that the platform crust is thicker than the shield, we find no statistically significant difference between whole crustal thickness of the Arabian shield and platform. However, the average sub-sedimentary crustal thickness (i.e., the crystalline crust) of stations on the platform is 34 km, 3 km thinner that the crust of the shield. Individual station Vp/Vs wavespeed ratios are highly variable for the Arabian plate, ranging from 1.60 to 1.97 and averaging 1.75, with a standard deviation of 0.07. There are no statistically significant differences between Vp/Vs ratios of the different geologic regions of Arabia. Similar Vp/Vs ratios, coupled with similar crustal thicknesses for harrats and shield, imply that Cenozoic magmatism has contributed negligibly to crustal growth.


Introduction
The present-day thickness and composition of the Earth's crust and its lateral variability are the product of its evolution from initial formation and evolution through multiple tectonic episodes. To infer the thickness of the present-day crust with minimal assumptions, we computed P-wave receiver functions (PRFs) and measured the travel-time differences between seismic P-waves and the converted phases (Ps, 3p1s, 2p2s) they generated at the boundary between the crust and mantle (i.e. the Moho). Compared to previous PRF studies of crustal thickness of Saudi Arabia (Sandvol and others, 1998;Julià and others, 2003;Al-Damegh and others, 2005;Tkalčić and others, 2006;others, 2016, 2019), we: (1) use more stations that were active for a longer recording period than in previous studies; (2) more thoroughly explore the methodological parameter space to obtain estimates of uncertainties; (3) correct for sedimentary cover where known; and (4) make statistical comparisons between different geologic regions.
The Arabian shield constitutes the eastern portion of the Arabian-Nubian shield (ANS) and is separated from the Nubian shield by the Red Sea rift, a divergent plate boundary that merges into the Dead Sea Transform boundary at its northern end ( fig. 1). The ANS was formed by the accretion of multiple terranes, primarily in the period~800-550 Ma (Stern and Johnson, 2010). The Arabian craton began to subside at 725 Ma and gradually accumulated up to 15 km of sediments (Stern and Johnson, 2010). The arrival of the Afar Plume at 30 Ma triggered extension in East Africa, the Gulf of Aden, and the Red Sea region (Stern and Johnson, 2010). Voluminous flood basalts in Ethiopia and Yemen were succeeded in Arabia by modest alkali-basalt fields, locally known as harrats (Stern and Johnson, 2010). Recent magmatism was recorded in the 2009 intrusion at Harrat Lunayyir and eruptions in Mohammedan times of Harrats Rahat and Khaybar, which lie along the Makkah-Medina-Nafud (MMN) line ( fig. 1a) (Camp and Roobol, 1992;Pallister et al., 2009).
In this paper we introduce the newly expanded seismic data set, then discuss the calculation of our receiver functions including data quality control. We describe the standard H-k stacking approach and our four key improvements over previous work in Saudi Arabia: (1) careful analysis of the effects of the chosen stacking weights; (2) use of a resonance removal filter (Yu and others, 2015), followed by a sedimentary H-k stack on stations whose data are difficult to analyze due to sedimentary reverberations; (3) application of phase-weighted stacking (PWS) (Crotwell, 2007) to noisy data; and most importantly, (4) application of an additional sedimentary correction based upon the well-mapped sedimentary thicknesses of the Arabian Platform (Stern and Johnson, 2010). We discuss our statistical approach in determining if apparent differences between subsets of the data are real, and then provide a geographic tour of our results.

Preprocessing and Quality Control
We examined data from 165 broadband seismic stations of the Saudi Geological Survey (SGS) National Seismic Network (Endo and others, 2007) that are now deployed in all of Saudi Arabia except for the Rub' al Khali (Empty Quarter) ( fig. 1a). In addition to the SGS data, we have also incorporated nine seismic stations from a 1995-1997 IRIS/PASSCAL array (Vernon, 1995). The recorders used by both SGS and IRIS are Trillium T40, T120, and Streckeisen STS2 seismometers (Blanchette and others, 2021). We obtained useful results at 154 stations (fig. 1b) from 898 teleseismic events ( fig. 1c) over a two-year period from 1995-1997 and a six-year period from 2008-2014. Each three-component recording was used to construct a P-wave receiver function (PRF) (Langston, 1979). The PRF method is a well-established technique to image seismic boundaries within the lithosphere, including the Moho and the base of sedimentary rocks that overlay the crystalline crust.
We rotated event recordings from their North-East-Vertical (NEZ) coordinate system to the Radial-Transverse-Vertical (RTZ) coordinate system, based on the theoretical back-azimuth of each recorded teleseismic earthquake, in order to isolate the first arrival (P-wave) on the vertical component and the converted S-wave on the radial component. We picked arrival times for every station-event pair using a short-term average (STA) over long-term average (LTA) algorithm (Earle and Shearer, 1994), and manually refined these arrival times to correct for slight errors in the picking algorithm by visually inspecting the data.
We band-pass filtered the raw data from 0.2 to 1.5 Hz and trimmed the filtered traces such that each trace begins 20 seconds before the first P-wave arrival and extends to 80 seconds after. PRFs were calculated by deconvolving the vertical (Z) signal from the radial (R) signal using an iterative-time domain approach (Ligorria and Ammon, 1999). Some previous workers simply visually inspected the calculated PRFs for quality (Sandvol and others, 1998;Thurner and others 2015;Miller and others, 2018;etc.); Tang and others (2016) used a minimum-fit criterion; and we quantitatively used both a minimum-fit criterion and a cross-correlation filter.

Figure 1. a) True-color image of Saudi Arabia
We band-pass filtered the raw data from 0.2 to 1.5 Hz and trimmed the filtered traces such that each trace begins 20 seconds before the first P-wave arrival and extends to 80 seconds after. PRFs were calculated by deconvolving the vertical (Z) signal from the radial (R) signal using an iterative-time domain approach (Ligorria and Ammon, 1999). Some previous workers simply visually inspected the calculated PRFs for quality (Sandvol and others, 1998;Thurner and others 2015;Miller and others, 2018;etc.); Tang and others (2016) used a minimum-fit criterion; and we quantitatively used both a minimum-fit criterion and a cross-correlation filter.
Our minimum-fit criterion requires that we are able to reproduce the radial and transverse components of motion with a fit that is equal to or better than 70% when we convolve the respective components of the PRF with the vertical component. The 70% minimum-fit criterion is chosen to enforce a decent match between the original signal and the PRF while not discarding an excessive amount of data. The initial dataset comprised 22,716 receiver functions, and 18,242 (~80%) of the receiver functions met our 70% fit requirement. Our crosscorrelation filter constructs a template PRF by stacking all PRFs for a given station and calculates the maximum cross-correlation coefficient of each event with the template. To minimize assumptions about the structure beneath each station, no moveout correction was applied to the PRFs, and we required a low minimumthreshold cross-correlation coefficient of 0.6 on the radial component of the PRFs only. Of the 18,242 receiver functions that passed the minimum-fit criteria, 15,519 (~68% of the raw dataset) passed the cross-correlation filter requirement. The final number of receiver functions is far lower on the platform (~11% of the final usable PRFs) than on the shield (~24%) or in the harrats (~51%; fig. 1b). A very high proportion of the useable receiver functions was generated by Mw 6+ subduction zone earthquakes northeast of the recording array (fig. 1c) in comparison to the less numerous and smaller mid-ocean ridge events at most other back-azimuths.

Standard H-k Stacking
PRFs utilize teleseismic P-wave energy from earthquakes that are 30-90°distant from the receiver. Some of the P-wave energy is converted to transmitted S-wave motion at first-order (abrupt) impedance boundaries in the lithosphere. This S-wave is denoted Ps (capital letters denote waves that are initially downgoing and lowercase letters indicate purely up-going waves). Three other arrivals ('multiples') are also commonly observed on the better-quality P-wave receiver functions: the PpPs (or 3p1s) arrival and two phases that arrive together (for a 1D Earth model) with the same sense of motion (PsPs and PpSs, collectively referred to as 2p2s) ( fig. 2). The piercing point at the Moho varies between phases and with the earthquake distance, but is farthest away from the seismic station (about 50 km, fig. 2c) for the 3p1s multiple at the shortest teleseismic offset used, here 30°. For clarity, we denote the depth of the interface at which the waves are refracting and converting with a subscript (e.g. P35s is a P-to-S conversion at 35 km depth). The travel-times (t [s]) of each phase as a function of crustal thickness (H [km]), P-wave velocity (Vp [km.s -1 ]), Vp/Vs ratio (k), and rayparameter (p [s.km -1 ]) are given by equations (1)-(3): Figure 2. Schematic ray paths of P-wave receiver function phases from a shallow earthquake at 30°d istance, for a 35-km-deep Moho separating crust (Vp = 6.1 km.s -1 , k = 1.73) and upper mantle (Vp = 8.04 kms -1 , k = 1.80). a) Direct P-wave. b) Ps phase. c) PpPs (or 3p1s) phase. d) PsPs (or 2p2s) phase. e) PpSs phase (also a 2p2s phase). Red lines are the labelled phase; black lines are the phase(s) in the previous panels.  H-k stacking (Zhu and Kanamori, 2000) is a standard approach used to analyze PRFs for H and k beneath a given seismic station, normally assuming a homogenous isotropic layer over a half-space beneath the seismic station. Standard H-k stacking employs a grid-search over H-k space for the stacking amplitude: where the rn are the amplitudes of the n th radial receiver function for this station at the predicted phase arrival times for the appropriate (H, k, Vp, p), and wi are the stacking weights for each phase. The normalization condition applied to the weights is 1 + 2 + 3 = 1, and we describe H-k stacks and their results as 'w1/w2/w3' to emphasize the weights used. The third term represents the 2p2s phases and is subtracted because it has the opposite polarity with respect to the other phases. The sum, s, reaches its maximum value for any seismic station at the average crustal thickness and velocity ratio beneath that station. We searched over a grid spanning 20 ≤ H ≤ 80 km and 1.60 ≤ k ≤ 2.10 to sample all plausible crustal thicknesses and compositions, using an assumed average crustal wavespeed (Vp avg = 6.5 km.s -1 ; Mooney and others, 1985). H-k estimates only vary weakly with Vp avg , increasing H by <~10 km per km.s -1 and k by <~0.05 per km.s -1 (e.g. Karplus and others, 2019). Conventionally the uncertainty in H and k values is computed using: where , , and ㄹ are the standard deviations of H, k, and the stack amplitude evaluated at the H and k values that maximize the stack (Zhu and Kanamori, 2000). We cite uncertainties of ±2 throughout this paper for individual stations, and in Tables, except in the Discussion section, where we cite uncertainties of ±1 for ensemble means.
It is common to use a default set of stacking weights, regardless of the signal-to-noise ratios of the different phases at different stations. Instead, we examined the individual contribution of each converted phase at each station to determine the final weights, selecting from among three options (0.4/0.3/0.3; 0.33/0.33/0.33; 0.5/0.5/0.0) to give the visually best focused result for the Moho conversions (cf. Ogden and others, 2019). Additional combinations were used to stack shallower conversions, within or at the base of the sedimentary cover, as discussed below. For example, at shield station AFFS ( fig. 3), each phase ( fig. 3c-e) provides clean constraints in H-k space. We generated H-k solutions for multiple combinations of weights and compared the predicted moveout curves of each phase for each H-k solution with the data to ensure the H-k solutions represent stacked coherent signal, not noise ( fig. 4). For station AFFS, our preferred weights are 0. 4/0.3/0.3 (fig. 3f).
In contrast to AFFS, at coastal station MWLHS ( fig. 5), the 2p2s phases provide potentially misleading information ( fig. 5e), and our preferred stacking weights are 0.5/0.5/0.0. We illustrate this issue in figure 6: the 2p2s stack (fig. 5e) suggests possible Moho depths of about 21 km or about 43 km. However, the predicted 2p211s and 2p431s arrivals (where subscripts denote depth of model converter) correspond to unusually high amplitudes on the more numerous PRFs with a small slowness of p~0.45 s.km -1 (the most distant earthquakes from the northeast) and not to any corresponding negative (blue) amplitudes on the sparse PRFs with larger p. Hence, we ignore the 2p1s stack in selecting our final weighting of 0.5/0.5/0.0 (Fig. 5f). It should be noted that although standard deviations calculated via equations (5) and (6) are appropriate for well-resolved H-k plots with sub-equal contributions from the various phases (e.g. figs. 3f and 5f), they are gross under-estimates for poor data quality particularly where only a single phase is available or dominates the stacks (e.g. figs 3c-e and 5c-e). For AFFS, our preferred result (H = 35.2 ± 0.4 km, k = 1.75 ± 0.02) is 8 H and 6 k from the value given by only the Ps phase ( fig. 3c).

Multiple removal to improve H-k Stacking
Standard H-k stacking works well at recovering the thickness of the crust, provided that the Moho is the only layer with a sharp enough impedance contrast to generate strong P-wave to S-wave conversions (Langston, 2011;Yu and others, 2015) and that the crust can be reasonably well-characterized with a single average wavespeed. However, some stations on the Arabian Platform are underlain by as much as 10-km of low-wavespeed sedimentary rocks ( fig. 1b) (Stern and Johnson, 2010), producing inter-bed conversions and multiples (reverberations) that can interfere with the Moho Ps arrival and its subsequent multiples. A significant number (21 out of 174) of the seismic stations analyzed here were affected by sedimentary multiples and required additional analysis.
We attempted direct removal of the sedimentary reverberations from these PRFs using a resonanceremoval filter calculated from the auto-correlations of the PRFs (Yu and others, 2015). This method estimates the apparent thickness and Vp/Vs ratio of the sedimentary layer, as well as a thickness and Vp/Vs ratio for the sub-sedimentary crystalline crust. This approach works best for sedimentary layers >~0.25 km thick for which the reverberations have frequencies and travel-times that interfere with the Moho conversions. However, this method may be inapplicable for layers >~8 km for which the wavespeed contrast at the sedimentary/basement contact is often sufficiently small or gradational that the reverberations do not have a high amplitude. We successfully filtered 15 stations with 2.6-10.7 km of sedimentary rock (thicknesses from Stern and Johnson, 2010). We then determined the crustal thicknesses using the H-k method on the filtered PRFs. As in conventional H-k stacking, to get a sedimentary thickness hs we assume a sedimentary Pwavespeed. We use 4.0 km.s -1 , based on refraction-profiling results showing Vp = 2.0-5.7 km.s -1 in the northern platform (Seber and others, 1993;Brew and others, 1997) and 3-5 km.s -1 within the western platform (Mooney, 1984), and carry out a grid search over H-k space. The sedimentary thickness, hs, is reasonably well determined by the Ps conversion from the sediment/basement interface, but determining the sedimentary Vp/Vs ratio, ks, requires good constraints on the corresponding 3p1s and 2p2s multiples that tend to have low amplitudes.

Phase-Weighted H-k Stacking
At five stations on the Red Sea rift margin (i.e. the coastal plain) for which neither the standard H-k stacking nor the resonance-removal filter proved useful, we successfully applied Phase-Weighted Stacking (PWS) (Crotwell, 2007;Ogden and others, 2019) to remove incoherent noise from the data. In this method, the coherency between the converted (Ps) and multiple (3p1s, 2p2s) arrivals is used to modulate the stacking amplitude s(H,k) : 9 where n are the instantaneous phases of the PRF (derived from the Hilbert transform of the PRF) at the predicted Ps, 3p1s, and 2p2s arrival times for the appropriate (H, k), and v is an arbitrary exponent. Using v = 0 corresponds to regular H-k stacking, with v = 1 (Crotwell, 2007, and this chapter), and v = 2 (Ogden and others, 2019) representing progressively more aggressive PWS. PWS uses the relative coherence of the arrivals to determine their contribution to the H-k stack. If the instantaneous phase  of the PRFs at the three predicted arrival times is equal (as expected for noise-free data), the weighted sum of the PRFs at these three times is given full weight in the H-k stack, whereas if  at the three arrival times is contradictory (as one would expect from random noise), the weighted sum of the PRFs is given much lower, or zero weight. In our application of PWS, we did not try to select optimum weights w for Ps, 3p1s, and 2p2s, but rather assigned each phase a value of 1/3. We attempted PWS H-k stacking only for stations for which we could not otherwise obtain reliable results. We obtained reliable results from an additional six stations. The remaining 14 stations defied our efforts to obtain a reliable crustal thickness, due to either the very low signal-to-noise ratios and/or very few PRFs at these stations.  Figure 4.

Sedimentary corrections to crustal thickness
Our single-layer interpretations of our H-k stacks (whether conventional or using PWS) assume a constant wavespeed for the entire crust. However, the traveltime of the Ps phase per kilometer of sedimentary rock (in a lower-Vp, higher-k basin) is larger than in the higher-Vp, lower-k crystalline crust. Hence, the presence of a basin makes crustal thicknesses directly measured from H-k stacking (Hraw) too large by an amount that increases in proportion to the basin thickness. For a sedimentary layer thickness, hs above a crystalline crust thickness hb (such that true crustal thickness Hcorrected = hs + hb), the travel time of the converted Ps phase with respect to direct P is given by: where the subscripts "b" and "s" denote crystalline basement and sedimentary values; so Hraw is given by: where subscript "c" denotes the assumed whole crustal values for the P wavespeed and the Vp/Vs ratio used in the H-k stacking. To correct Hraw to the true Moho depth H, we use the sedimentary thickness hs (equivalent to the basement depth) that is well-known across the Arabian platform from extensive hydrocarbon exploration (e.g. Stern and Johnson, 2010) ( fig. 1b). We define the difference between the true crustal thickness and the raw thickness as the error, f = 㐰 − . For PRFs the ray-parameter, p, ranges from 0.04 s.km -1 to 0.08 s.km -1 , and in figure 7a, we show the variation of f with sedimentary thickness hs (i.e. f ℎ ㄹ ) as a function of sedimentary Vp and Vp/Vs ratios for constant p = 0.06 s.km -1 . For a sedimentary layer thickness of hs, with Vp = 4.0 km.s -1 and ks = 1.75 (Brocher, 2005; fig. 7), Hraw must be reduced by~0.6 hs. Figure 7b shows the variation of f with slowness p (i.e. f ) is irrelevant for the range of p used here. In the following sections, we present the full crustal thickness H (Moho depth below the surface), after reduction by 0.6 hs using sedimentary thicknesses from Stern and Johnson (2010), and we also discuss crystalline crustal thickness hb = H -hs. At stations for which we performed the resonance-removal and sedimentary stacking method (Yu and others, 2015), we applied the 0.6 hs reduction to the difference between the H-k sedimentary thicknesses and the known sedimentary thicknesses. This correction was applied because, in general, the H-k sedimentary thicknesses are less than the known sedimentary thicknesses, meaning we are still over-estimating the total crustal thickness.
There are exactly two stations for which we obtained sedimentary H-k values greater than the thicknesses we expected, therefore, at those two stations, the correction is in the opposite direction.

Comparison of H-k to common-conversion-point (CCP) processing
The importance of the sedimentary correction is shown by comparing H-k results before and after correction, overlain on common-conversion-point (CCP) images (fig. 8) along Red Sea-"parallel" and "orthogonal" profiles ( fig. 1a). The CCP stacks were calculated on a latitude/longitude grid with 50-km spacing between nodes in each horizontal direction and 2-km spacing vertically. Each CCP node had a 50-km radius, over which rays were stacked, i.e. rays parallel to each profile are used by two adjacent CCP nodes; orthogonal rays are used only once. Hence the alignment of the "orthogonal" profile towards the greatest earthquake density ( fig. 1c) is one reason the "orthogonal" CCP image ( fig. 8b) is smoother than the "parallel" CCP image ( fig. 8a). To create grid-oblique CCP cross sections, we chose only grid nodes within 25 km of the transect. Points along the transects that had more than one node within the 25-km distance limit are composed of the average of the nodes. On the Red Sea-parallel section, the H-k stacking and CCP stacking results are mostly in good agreement, including in places where there are very large abrupt shifts in imaged Moho depth, e.g. near 24°N , where the station layout along our linear profile switches from stations only on the shield (Moho at~40 km) to a station only on the coastal plain (~20 km Moho) and back on the shield again ( fig. 8a). Minor discrepancies likely result from assuming the same 1D velocity model across all stations for CCP stacking, whereas H-k stacking solves for a unique layer over a half-space velocity model at each station. An additional source of discrepancy comes from CCP stacking using multiple stations at some nodes, whereas H-k stacking is purely single-station. Large differences (e.g.~27.5°N) appear due the inability of CCP stacking to handle multiple candidate Ps peaks, whereas H-k stacking enforces consistency between the Ps, 3p1s, and 2p2s phases. The Red Sea-orthogonal section shows a consistent sense of offset, with the H-k Moho shallower than the CCP stacking Moho by up to 8 km where the sedimentary section is thickest along our profile, as expected from our analysis as described above. In spite of the apparently well-resolved Moho on the CCP image, it is clear that care must be taken in interpretation of such an image due to the additional travel-time delay of the sediments.

Results
We next show H-k analyses corroborated by direct inspection of P-wave receiver functions, including examples from and compiling crustal-thickness results for each of the geological regions (Harrats, Coastal Plain, Shield, Platform; fig. 1a).

Harrat Lunayyir
Harrat Lunayyir ( fig. 1a, 9a) has the greatest station density of the SGS network, with 18 seismometers within an~50-by-50 km area (Blanchette and others, 2018; this volume), which were installed during and following the 2009 dike intrusion. This dike intrusion was accompanied by >30,000 earthquakes, up to Mw 5.7 (Pallister and others, 2010). One station (LNY07) in this sub-network is <250 km from the Red Sea rift axis and has elevation < 400 m above sea-level and is here considered to be a 'coastal' station.
LNY02 is a representative seismic station near the center of Harrat Lunayyir ( fig. 9a). The PRFs at station LNY02 ( fig. 10a) show a coherent arrival that we identify as Ps, with a positive amplitude (red), arriving 4-5 s after the initial P-wave arrival. Assuming a simplified 1D velocity model (Vp = 6.5 k.ms -1 , k = 1.80), with zero sediment thickness (hs) and 34-km Moho depth, we can reasonably match the travel-times of this arrival as P34s (Moho conversion from 34 km depth), and the later positive arrival at a delay-time of~14 s as P34pP34s, i.e. the 3p1s Moho multiple. In practice, we examined the H-k stacks ( 9d) provides a different but still well-defined zone of possible solutions in H-k space, that is also insensitive to H and, particularly, k. Although the 2p2s multiples can be used to form a H-k stack ( fig. 9e), this does not conform to the other two phases, and is seen to correspond to a diffuse region of negative and low-amplitude energy on the PRF ( fig. 10a). We chose stacking weights 0.5/0.5/0.0 at this station and obtain H and k values of 34.0 ± 0.6 km and 1.80 ± 0.04. These values match those of Tang and others (2016), who reported H = 34.1 ± 0.4 km and k = 1.80 ± 0.02 (2 uncertainties). Results from equivalent analyses of all Harrat Lunayyir stations are given in Table 1, with average Moho depths of 33.7 km and crustal k of 1.78. In a 1D and isotropic Earth model, we expect the transverse component of the PRF to be composed purely of noise (Savage, 1998;Schulte-Pelkum and Mahan, 2014), so coherent transverse energy ( fig. 10b and c) likely represents an obliquely dipping boundary beneath the station. It is possible that such a dipping structure and corresponding lateral velocity heterogeneity has defocused the 2p2s multiple ( fig. 10a). Such considerations suggest a heuristic uncertainty of at least ±1 km in Moho depth and ±0.05 in Vp/Vs ratios. The coherent arrival between the Ps and 3p1s phases (6-7 s) cannot be explained as a conversion or a multiple from the Moho; nor is it a sidelobe of the Moho P34s, behind which it is delayed by ≥ 2s, clearly more than the central period of our data (T/2 < 1.5 s). The negative amplitude of this arrival is, however, consistent with a P-to-S conversion from a deeper, negative, velocity discontinuity with wavespeed decreasing downwards. We identify this arrival as the Ps conversion at the lithosphere-asthenosphere boundary LAB (Blanchette and others, 2018;this volume). The incoherence of the expected 2p342s phases is not due to superposition with the LAB-generated phases ( fig. 10a). The seismic velocity contrast across the LAB is thought to be of smaller magnitude and spread over a wider vertical transition zone than the contrast across the Moho (Artemieva, 2011). These factors would yield a lower-frequency, lower-amplitude signal than the Moho conversion, the reason, we believe, we did not observe multiple arrivals from the LAB on this or any other station. Although this LAB arrival is evident on numerous stations, in this paper we reserve attention for the Moho.
There is one Harrat Lunayyir station, LNY06, for which H-k stacking is ambiguous. LNY06 is only~12 km from station LNY02 ( fig. 11a), but has very 'ringy' data (many coherent phases, possibly multiple reflections from within or at the base of basalt flows) ( fig. 12a), and shows two distinct peaks in H-k space (fig. 11f) at depths of 30.4 ± 0.5 km and 34.4 ± 0.5 km. Although we would be unable to pick the correct Moho conversion from this station alone, the other 17 nearby stations over Harrat Lunayyir have 31.2 < H < 37.4 km, so that the most plausible result at station LNY02 for H is 34.4 ± 0.5 km, which we report here.  LNY06 and LNY02 both appear to show a converter from the LAB but at delay times respectively of 6.25 s and~6.5 s, suggesting variable LAB depth (figs. 10, 12).  Figure 3.

Harrat Khaybar
Harrat Khaybar is the northernmost harrat along the Makkah-Medina-Nafud (MMN) line (Camp and Roobol, 1992) (fig. 1) of alkalic volcanism that was emplaced since~12 Ma during the second stage of Red Sea rifting (Blanchette and others, 2018). Whereas prior to 2012 there was only one seismometer on the western edge of Harrat Khaybar (KBRS), there is now an array of 10 stations within Harrat Khaybar ( fig. 1).
The seismometer KBR08, on the western edge of Harrat Khaybar ( fig. 13a), has PRFs ( fig. 14a) that share many similarities with the observations discussed for Harrat Lunayyir stations ( fig. 12a). There is a coherent, positive-amplitude (red) Moho Ps phase at~4-5 s delay time, followed by a negative-polarity LAB Ps phase at~7 s. The 3p1s phase is coherent and arrives between~14-15 s. As for LNY02 there is a hint of 2p2s at~17-20 s, but the 3-s time span and low amplitude suggest lateral wavespeed heterogeneity. Significant Moho dip is ruled out by the consistent phase of Ps at all back-azimuths (not shown; see Blanchette and others, 2021) (Schulte-Pelkum and Mahan, 2014).
These observations from the KBR08 receiver functions are consistent with its single-phase H-k stacks ( fig. 13 c-e). The Ps stack ( fig. 13c) and 3p1s stack (fig. 13d) are well resolved, but the 2p2s arrivals ( fig. 13e) are too diffuse, spanning too broad a range of H to provide much sensitivity. We use 0.5/0.5/0.0 to obtain H = 36.4 ± 0.7 km and k =1.72 ± 0.02, identical within uncertainty to previous results (35.1 ± 0.6 km and 1.74 ± 0.04) for station KBRS, which is just 23 km distant (Tang and others, 2016). Our own analysis for station KBRS results in H and k values of 38.2 ± 0.5 km and 1.76 ± 0.02, emphasizing the distinction between formal and heuristic uncertainties. Results for all Harrat Khaybar stations are given in Table 2, with an average Moho depth of 35.5 km and crustal k of 1.75.

Harrat Rahat
Northern Harrat Rahat last erupted in 1256 C.E. (Camp and others, 1987) and was the focus of the 2013-2017 cooperative SGS-USGS project to investigate volcano and seismic hazards to the City of Medina. There are 14 seismic stations within the densely instrumented northern portion of Harrat Rahat ( fig. 1 and 15b). Three stations within this sub-network, RHT05, RHT06, and RHT13, are > 30 km (~one crustal thickness) from Cenozoic basalt flows (Pollastro, 1998a) and are here considered 'shield' stations; RHT16 and six other stations were installed in the southern portion of Harrat Rahat, the southern extent of the MMN line ( fig. 1a), and are included here as "harrat" stations. The receiver functions at this station ( fig. 16a) are as clear as those at KBR08. There is a broad, low-amplitude positive arrival between 0-2 seconds, with a moveout opposite to that expected from a shallow Ps conversion. This phase appears to have slightly smaller differential travel-times from events originating with a back-azimuth of~70-130°, unlike the Ps phase. The Moho Ps arrival and the 3p1s phase are both narrow and high-amplitude. The 2p2s arrivals are strongest at ray-parameters (p) ≤ 0.055 s.km -1 , likely due to our uneven earthquake source distribution (figs. 1b, 16a). These phases are consistent with the single-phase stacks ( fig. 15 c-e), and the stacking results are relatively insensitive to our choice of weights ( fig. 15f). We chose to use 0.4/0.3/0.3 as our stacking weights with resulting H and k of 35.0 ± 0.6 km and 1.75 ± 0.03, essentially identical to those of Tang and others (2016) who obtained H of 35.1 ± 0.4 km and a Vp/Vs ratio of 1.74 ± 0.02. Station RHT14 is located on the youngest lava flows close to Medina (fig. 17a). The P-wave receiver functions for this station are dominated by strong ringing that was not significantly reduced by the resonance filter ( fig. 18a) but still contain useful information about the crust. There is a positive arrival of unknown origin following the direct P-wave by < 1 s: the 1256 C.E. lavas are ≤ 100 m thick (Downs and others, 2018), so would produce a signal delayed by only~0.02 s, but these lavas obscure the deeper geology.
Like the direct P-wave, the Moho Ps (4-5 s) is also followed by complicated reverberations. There are weak and diffuse 3p1s and 2p2s arrivals that are better identified by their agreement with predicted arrival times than by their absolute amplitudes (fig. 18a). These 3p1s and 2p2s arrivals provide very streaky images in the single-phase H-k stacks ( fig. 17c, d). Although the analytical uncertainty would be lower if we excluded the 2p2s arrivals from the H-k analysis, we chose 0.4/0.3/0.3 stacking weights to incorporate the information in the 2p2s phase and obtained H and k values of 35.8 ± 1.0 km and 1.76 ± 0.03, respectively. The H-k stacking results for the rest of the Harrat Rahat stations are listed in Table 3. The average Moho depth is 35.2 km and crustal k is 1.75.

Red Sea Rift
A total of 29 SGS stations are within 250 km of the Red Sea rift axis and at elevations ≤ 400 m a.s.l. (fig.  1a). The Great Escarpment, a prominent topographic feature defining the Red Sea rift shoulder, is typicallỹ 100 km from the coast, and rises above 1,000 m in southern Saudi Arabia, becoming lower north of Jeddah and eventually indiscernible to the north ( fig. 1b). We include these 29 stations as lying within the Red Sea rift, including two stations close to the Gulf of Aqaba/Dead Sea Transform, and two stations on Farasan Island. The two northernmost coastal stations lie on the Dead Sea Transform (DST), the~1,000-km-long rightlateral strike-slip fault that stretches from the northern tip of the Red Sea to the East Anatolian Fault. The crust near the DST has been previously measured to be 32-37 km thick from wide-angle reflection/refraction experiments (El-Isa and others, 1987;DESERT Group, 2004;ten Brink and others, 2006;Mechie and others, 2009) and PRF analyses others, 2005 and. Because of the extensive previous work on this structure, we simply report results from the two SGS seismic stations~15 and 50 km from the DST, HAQS, and BDAS (Tables 4 and 5). We show HAQS (Figs. 19, 20) as exemplary of these northernmost stations, even though we list it in Table 5, as a Shield station, due to its elevation of 450 m. At station HAQS, we chose optimal weights of 0.4/0.3/0.3 and obtain H and k values of 28.0 ± 0.7 km and 1.77 ± 0.03, respectively. Station QNF01 lies on the coast of the southern Red Sea rift margin ( fig. 21a) and is one of the stations at which we were unable to obtain a reliable estimate of H and k via standard H-k stacking, as both the H-k plots ( fig. 21c-e) and the PRFs ( fig. 22) are noisy. Instead, we used PWS H-k stacking ( fig. 21f) that requires consistency between the Ps ( fig. 21a), 3p1s (fig. 21b), and 2p2s (fig. 21c) phases and yielded H = 36.6 ± 0.3 km and k = 1.71 ± 0.01. Moveout curves for each of these phases are plotted in figure 22a and appear to match coherent signal in the PRFs, even though the resulting crustal thickness of 36.6 km indicates a surprising lack of crustal thinning here.
Tertiary sedimentary rocks are locally more than 4 km thick at the coast (e.g. Cole and others 1995), but most of the seismic stations are deployed on or close to crystalline basement and were not corrected for overlying sediments. The probable exceptions are QNF01 ( fig. 21 and 22) in the Ghawwas sub-basin (for which phase-weighted stacking was required to achieve a solution), which we corrected for an assumed 1-km sedimentary thickness (Hughes and Johnson, 2005), and the Farasan Island stations (FRSS and FRSS2), which we corrected for 0.6 km of sedimentary rocks (Almalki and Bantan, 2016). Table 4 gives H-k results for all the 'rift' or 'coastal' stations.  Table 5.

Platform
Platform stations have predictably worse data quality than shield stations (due to higher noise levels compared to stations on hard rock), with only 12% of shield stations (but~25% of platform stations) not producing an acceptable H-k result ( fig. 1b). Data quality typically worsens, and the number of acceptable PRFs decreases, as the sedimentary thickness increases to the east and intra-sedimentary reverberations become stronger ( fig. 1b) 5 s (fig. 24) is from a very shallow converter, far shallower than the~7-km Phanerozoic sedimentary thickness known in this area ( fig. 1b). This example illustrates the need to correct our raw crustal thicknesses for the entire Phanerozoic basin thicknesses taken from Stern and Johnson (2010) (fig.  1b). H-k stacking results for all platform stations are given in Table 6 below.

Regional Averages
Our H-k stacking results ( fig. 25a & b) show that crustal thickness of the Arabian plate is least near the coast of the Red Sea and increases away from the rifted margin, as expected. The average total thickness of the crust across the entire array is 34.9 ± 4.7 km (1 uncertainty). The average total thicknesses of the rift-margin (coastal), shield, and platform crust (29.1 ± 4.1 km, 36.6 ± 4.8 km, and 38.4 ± 3.4 km, respectively) are all close to the global averages for these crustal types (30.5 ± 6.0 km, 40.8 ± 7.0 km, and 39.0 ± 7.0 km, respectively) (Mooney and others, 1998). Analytic errors underestimate the uncertainty in the value for the Arabian platform because of the inability of the sedimentary H-k stacking to obtain true sedimentary thicknesses and the necessity to assume an a priori sedimentary Vp when making corrections for sedimentary thickness (as discussed in our methods section). The average total crustal thickness of the harrats is 34.7 ± 1.6 km.           Corrected thickness = H -0.6 • hs-SJ2010 where hs-SJ2010 is sedimentary thickness from Stern and Johnson (2010) plus elevation a.s.l.; where hs is available from sedimentary H-k stacks, corrected thickness = H + hs -0.6 • (hs-SJ2010 -hs).
In map view, the crustal thicknesses show systematic variation ( fig. 25 a, b, c) that is not evident for the Vp/Vs ratios ( fig. 25d). The average Vp/Vs for the entire array is 1.76 ± 0.07, deviating only slightly from a Poisson solid (1.73). Averages for each tectonic province are given in Table 7. Although crustal thicknesses have distinct means ( fig. 26a), all regions of the Arabian plate have statistically equivalent k-values ( fig. 26b).  in the same provincial framework as we have done here. We reorganized their results (Table 8) following our classification scheme to enable comparison in sequence: coastal (Red Sea rift flank), harrat, shield, and platform. For coastal stations, Tang and others (2016) report average total crustal thickness of 27.9 ± 4.0 km across 11 stations and Al-Damegh and others (2005) 25.1 ± 6.1 (5 stations), consistent with our average of 29.1 ± 4.1 km. Tkalčić and others (2006) reported 28.0 km for a single coastal station (YNBS), for which we obtained 24.0 ± 0.7 km. This significant discrepancy is due to the greater focus of Tkalčić and others (2006) on fitting the long-period P-wave receiver function, instead of the short-period P-wave receiver function (as can be seen in their Figure 11 part j); for this station Tang and others (2016) reported 25.0 ± 0.8 km, close to our own result.
All of our average harrat results are consistent with the previous studies, likely because of large station numbers, generally high data quality, and no sedimentary interference. Our shield result (crustal thickness of 36.6 km) is lower than most, but within error of all, previous studies, likely due to our inclusion of a greater number of seismic stations in the western shield, closer to the coast. On the Arabian platform our average raw total crustal thicknesses (40.4 ± 4.3 km) is within error of previous results. Inclusion of the sedimentary correction (our total crustal thickness, average 38.4 ± 3.4 km) increases our discrepancy from the conventional analysis of Tang and others (2016) (42.4 ± 5.6 km), but causes convergence with Tkalčić and others (2006) (36.0 ± 2.6 km), suggesting that their inclusion of long-period surface waves in a joint inversion with PRFs, reasonably compensated for the sedimentary cover, without necessitating an additional sedimentary correction. No previous authors published estimates of crystalline crustal thickness (our average is 33.8 ± 3.6 km).

Statistical comparisons of regions
A long-standing question about the Arabian Plate is whether different crustal terranes have different crustal structures. Stern & Johnson (2010) claim that the platform is 3-4 km thicker than the shield, and Tang and others (2016) claim they can distinguish between crustal thicknesses of different Proterozoic terranes and that k-values at Harrat Lunayyir are anomalously high. However, these claims were not evaluated statistically. Here, we use different statistical methods to test the equivalence between different geographic regions.
First, we use Welch's t-test method (Welch, 1947), which generalizes the Student's t-test to allow the measurements being compared to contain both different sized populations and different variances. We use a two-tailed test of two scenarios (H0 -our null hypothesis, that the population means are equal, or H1 -the population means are different) and calculate the test statistic, t, and the degrees of freedom, , associated with the estimated variances. Calculating the Student's distribution with these input values allows us to calculate the probability, p, that the null hypothesis (H0) is true. At a 95% confidence threshold, if p < 0.05 we reject the null hypothesis (i.e. the averages are too different to be explained by random chance) and if p > 0.05 we accept the null hypothesis (i.e. any differences between the averages can be explained by random chance).
As an additional test of equivalence between geographic regions we also use the Mann-Whitney U (MWU) test (Mann and Whitney, 1947), designed to compare population distributions instead of population means, and to cope with non-normal, e.g. skewed, distributions (Fay and Proschan, 2010). We use the same two-tailed test as for our t-test and calculate p, the probability that the null hypothesis is true. We denote the two different p-values as either pt (t-test) or pmwu (MWU test). showing the probability density function (pdf) of the Student distribution for the two datasets. Blue vertical bars are 95% confidence bounds; the red bar is the t-value of this comparison. If the red bar is within the blue bars, then the null hypothesis (H0) is accepted, otherwise H1 -that the distributions are different -is accepted. In this case the two population averages and distributions are very different (pt < 10 -4 , pmwu < 10 -4 ).
The thinnest crust, 16-25 km, is beneath the Farasan Island stations. The two measurements, separated by only 37 km, differ greatly (Table 4), perhaps because of complex overlying sedimentary structure (Almalki and Bantan, 2016), but bracket the previous active-source measurement (17.5 km; Mooney and others, 1985) and represent the transition from thin oceanic to thick continental crust. The next thinnest crust is beneath the stations along the Red Sea rift margin. The average thickness of crust beneath the coastal stations is 29.1 ± 4.1 km, or 29.8 ± 3.2 excluding the Farasan Island results. Even excluding the Farasan stations, this value is statistically different from the average thickness of 36.2 ± 3.8 km for the rest of the Arabian plate ( fig. 27). We note that our definition of "coastal" is necessarily arbitrary, even though geologically informed (within 250 km of the Red Sea rift axis and at elevations ≤ 400 m a.s.l.), and other definitions would lead to different numerical values. Crustal Vp/Vs ratio for each station. White circles are approximately located at the center of each harrat (Lunayyir, Khaybar, and northern Rahat) and have a radius of 25 km (the averaging radius for fig. 25, and the distance-scale over which we anticipate commonality between stations due to similar ray path coverage ( fig. 2)). Thin black lines: harrat boundaries. Yellow star is Medina.
Moving eastward into the harrats, crustal thickness of Harrat Lunayyir is 33.8 ± 1.8 km and the Vp/Vs ratio (k) is 1.77 ± 0.09 (1). The scatter in k values far exceeds the analytic uncertainty of individual stations; however, this scatter exists between stations that are sufficiently close ( fig. 28) that they share many ray-paths contributing to their H-k results ( fig. 2), and should have similar results. Hence, we conclude that the low analytic uncertainties fail to represent the oblique eccentricity of resolution in our H-k grid search, which we represent as 90% amplitude contours. Because crustal H and k represent secular evolutionary processes, we might expect correlations within provinces ( fig. 29), as well as differences between them ( fig. 26b, 29a). Figure  29 demonstrates the lack of correlation between the thickness of the crust and its Vp/Vs ratio, except apparently amongst the harrat stations. However, the slope of the apparent trend of H versus k is the same as the slope of the H-k confidence ellipses ( fig. 29b), suggesting that the trend is due to real uncertainty in the H-k results, not geological differences. Indeed, adding mafic intrusions to increase k would also increase H. This is important because Tang and others (2016) proposed that the k-values at Harrat Lunayyir are elevated in comparison to the shield due to solidified intrusions. Our statistical testing methodology shows that the claim of Tang and others (2016) is not supported by their data, which consist of 11 Harrat Lunayyir measurements of k with a mean of 1.78 ± 0.11 and 36 shield measurements with a mean of 1.74 ± 0.07. Our t-test and MWU test applied to their data yields a pt-value of 0.35 and pmwu = 0.27, both far above the 0.05 threshold, so requiring the null-hypothesis that the Vp/Vs ratios of Harrat Lunayyir are statistically the same as for the rest of the shield.
Harrat Khaybar has average crustal thickness of 35.5 ± 1.3 km and average k of 1.75 ± 0.03. The consistency of results arises from both the low-noise levels of their receiver functions and the well-resolved phase arrivals. Harrat Rahat has mean H and k-values of 35.2 ± 1.1 km and 1.75 ± 0.05, respectively. The highest k-values for Harrat Rahat (≥1.78) are at the northern tip, where volcanics erupted~800 years ago (Camp and others, 1987), but, as discussed for Harrat Lunayyir, we believe this variability in Vp/Vs ratios is likely due to inherent uncertainly. The three major harrats in this study (Lunayyir, Rahat, and Khaybar) all have overlapping H and kvalues ( fig. 30). The k-values for the three harrats are not statistically distinct ( fig. 30c), but Harrat Lunayyir is distinctly thinner than Harrats Khaybar and Rahat, which are essentially identical ( fig. 30b). We suspect this thickness difference can be explained by the Harrat Lunayyir stations being closer to the Red Sea rift axis than all Khaybar stations and the majority of Rahat stations. We expand upon this overall west-to-east crustal thickening trend below, but the trend can be seen by comparing harrat versus shield crustal thickness (statistically distinct, fig. 26a), particularly with MMN versus western shield stations (statistically the same, fig.  31b).
Next, we compare the western shield (Htotal = 36.0 ± 5.4 km) to the MMN line (harrats Rahat and Khaybar; Htotal = 35.3 ± 1.2 km) ( fig. 31b). The western shield and the MMN line have the same mean thicknesses (pt = 0.44), but significantly different distributions (different variance) (pmwu = 0.04) which we reconcile by looking at the spatial distribution of measurements ( fig. 31a). The MMN line samples crust in a small region, so the measurements are tightly grouped around their mean. The western shield spans a much larger spatial region, including north to the DST and south nearly to the Yemen flood basalts, and thus have a more variable range of observations. This is the only example we found where the t-test and MWU test give different results, and is simply due to the different definitions of equivalence (same mean versus same distribution of values).   . 31c). The platform including its sedimentary cover may be thicker, but it is not statistically different (pt = 0.07, pmwu = 0.21). However, when we compare the thickness of the crystalline crust alone ( fig. 31d), the platform average is 33.8 ± 3.6 km, which is 2.8 km thinner than and distinct (pt << 0.05, pmwu << 0.05) from the shield. Although this result appears conclusive, it is heavily dependent upon how we correct for the sedimentary cover. We tested varying our sedimentary correction factor f ℎ ㄹ from -0.3 to -1.0 (corresponding to assumed sedimentary-basin average P wavespeeds of 3.5-4.8 km.s -1 ) ( fig. 7) and found that within this range, the platform crystalline crust is always statistically thinner than the shield (0 < pt < 0.05, 0 < pmwu < 0.01) ( fig. 31d). In the unlikely case that the sedimentary-basin average P wavespeed is greater than 4.8 km.s -1 the platform crystalline crust thickness becomes statistically indistinguishable from the shield.
Another potentially confounding factor is that the relative thicknesses of platform crystalline crust and shield crust may be biased by our station distribution. When we plot crustal thickness against distance from the Red Sea rift axis (RSR) ( fig. 32a), we see at the broadest scale that the Arabian Plate thickens from the western rifted margin to the eastern platform. The group of shield stations with H >40 km at distances of 200-400 km from the RSR does not fit this west-to-east thickening trend, in large part due to a tendency for our measured shield thicknesses to increase to the south ( fig. 25a, 32c). Total crustal thicknesses of more than 45 km close to the Yemen border may reflect the influence of plume magmatism that produced the Yemen flood basalts (Stern and Johnson, 2010). This increase in thickness to the south is also seen in the results of active-source seismic studies (Mooney and others, 1985) (fig. 1a) that show greater shield thicknesses than our averages ( fig. 32b). If we limit our comparison to only our southern stations that are within 125 km of the refraction profile ( fig. 32b), this discrepancy largely disappears. A best-fit straight line to the crustal thickness of all stations greater than 250 km from the RSR ( fig. 32c) shows the crust thinning northwards at 0.005 km/km. If this northward-thinning tendency is representative of the whole Arabian plate, then our bias in station distribution, with the centroid of our shield stations located~200 km south of the centroid of our platform measurements, would explain 1 km of the observed thickness difference between the crystalline crust of the platform and shield. It is possible that when more data become available for the Arabian Platform in southern Saudi Arabia (the Empty Quarter), the Arabian Shield and the Arabian Platform will be found to be statistically indistinguishable. We note that performing the same analysis on areal averages of the data (instead of station averages) does change the amount of weight different stations contribute to the analysis but does not change our final results. Our observation that crustal thickness increases with distance from the Red Sea Rift ( fig. 32a) remains apparent after correcting all data for 0.005 km/km northward-thinning ( fig. 32d). These data have the visual appearance of and can be fit with, an exponential curve ( fig. 32d). An exponential increase in thickness away from the RSR is appropriate for a plate cooling model, but we have no expectation that crustal thickness should be controlled by a rifting process far into the plate interior. We therefore also fit the data with two straight lines: a best-fit line of arbitrary slope (0.08 km/km) for all stations less than 250 km from the RSR; and a best-fit horizontal line (average thickness) for all stations greater than 250 km from the RSR. That these lines intersect at 250 km from the RSR validates our choice of this distance as the boundary between crust that has been thinned by Red Sea-rifting processes and the plate interior that is unaffected by rifting. The Great Escarpment ( fig. 1b) that is the geomorphic expression of the rift boundary is approximately 250 km from the RSR ( fig. 1a), suggesting that Red Sea rifting of the Arabian plate margin occurred via uniform crustal extension that affects a region no broader at depth than at the surface. Figure 32. Crustal thickness versus distance from and along the Red Sea rift (RSR) (defined with respect to axes in Figure 1a), displayed with crustal thickness increasing downwards, so plots are comparable with crustal cross-sections (Fig. 33). a) Crustal thickness versus distance from RSR. Triangles are singlestation crustal thickness estimates color-coded by region. For platform stations, we plot both total crustal (gray) and crystalline crustal (white) thicknesses. Squares are regional averages with 2 error bars. b) Same as part a) with crustal thickness values from Mooney and others (1985) (white stars), all stations > 125-km distant from Mooney et al. profile are now gray, and platform crystalline crust thicknesses no longer included. c) Crustal thickness versus distance along RSR, with distance increasing northwards.
Receiver functions contain valuable information on crustal properties, and extracting that information is straightforward beneath stations with low noise on crust that is relatively simple. We have extended this method to stations with higher noise and locations with thick surficial sediments by employing three steps: (1) examine each receiver function for quality and retain only the reliable ones, (2) evaluate the stability of the solution depending on different stacking weights for three well-recorded P-to-S-wave converted phases, and (3) employ more advanced forms of the H-k stacking algorithm. These three steps must be done in a self-consistent manner that honors the original waveform data. Finally, it is vital to correct crustal thickness estimates for the effects of lower-wavespeed sedimentary basins.