Constraining the Lithospheric Discontinuity Structure Beneath Hawaiʻi Using Teleseismic Receiver Functions

The Hawaiian Islands represent the youngest portion of the long‐lived Hawaiian‐Emperor Seamount chain sourced from hotspot volcanism on the northwesterly moving Pacific plate. Despite being one of the best studied hotspots on Earth, our understanding of the lithospheric discontinuity structure of Hawaiʻi is inhibited by a lack of a long‐term distributed broadband seismic network and lithospheric‐scale imaging of other oceanic hotspots. In this study, we analyze 3,530 teleseismic waveforms from 1,130 events recorded at 47 stations to compute P‐wave radial receiver functions. We then incorporate an elevation correction to adaptive common conversion point stacking to create a 3D discontinuity volume of the lithospheric structure associated with an ocean island hotspot. We find a positive amplitude conversion consistently at ∼12 km below sea level, likely representing the crust‐mantle transition below much of the island. However, this conversion occurs at only ∼6 km below Kīlauea, where local tomography images a zone of high P‐wave velocities (>7 km/s) within the crust and deformation of olivine phenocrysts within a magmatic reservoir is thought to occur. At this location, seismicity also extends vertically to 35 km depth before merging with an interpreted sill complex at ∼40 km. We interpret this anomalously shallow high amplitude discontinuity as an olivine‐dominated magmatic reservoir with an estimated maximum 5% melt over the resolution of our data at crustal depths. This represents the shallow seismic signature of a much more extensive magma system extending to depth. These images provide a better understanding of the magmatic structure associated with oceanic hotspot volcanoes.


10.1029/2023JB027029
2 of 13 The upper crustal structure beneath the island of Hawaiʻi has been investigated via geochemistry, seismicity, seismic imaging, and other geophysical methods making it relatively well understood compared to deeper structure.Variations in the geochemistry of erupted basalts are related to the maturity (age) of an oceanic hotspot volcano.Their corresponding magma supply rates and presence of xenoliths have been used to constrain the depth of magma reservoirs beneath active edifices from one to a few kilometers in depth (Bohrson, 2007).A synthesis of the separate magmatic systems underlying Kīlauea and Mauna Loa, the most active edifices today, uses geodetic data from multiple eruptions to identify an interconnected system with multiple reservoirs (Poland et al., 2014 and references therein) and agrees well with InSAR (Baker & Amelung, 2012).For example, multiple reservoirs of a few cubic kilometers in volume beneath the Kīlauea edifice are located at ∼1 km and between 3 and 4 km with a smaller magma chamber ∼1 km below the Eastern Rift Zone (ERZ).GPS measurements and waveform cross correlations indicate a degree of independence in the shallow (∼2 km) plumbing system between Kīlauea and its rifts; however, they appear to be fed by the same deeper reservoir (3-4 km) (Lin & Shearer, 2021;Neal et al., 2018).Mauna Loa, a more mature but similarly active edifice, shows some correlation with Kīlauea's eruption cycle, but the interconnectedness of their magmatic systems remains unknown as this dependence may be explained by stress transfer from the injection of new magma and subsequent eruptions (Przeor et al., 2022).
Other studies have used local passive and/or active seismic techniques to understand the relationship between intracrustal structure (<20 km BSL) and the crustal magmatic plumbing system beneath Hawaiʻi (Ellsworth & Koyanagi, 1977;Haslinger et al., 2001;Hill, 1969;Lin et al., 2014;Okubo et al., 1997;Park et al., 2009;Syracuse et al., 2010;Thurber, 1987;Zucca et al., 1982).These studies often feature high P-wave velocity regions (Vp ∼7 km/s) below edifices at 6-8 km depth that extend downward to the base of the crust (generally ∼15-20 km) or the base of the model depending on the depth sensitivity of the particular study.These high velocity regions are often interpreted as being mafic to ultramafic cumulates (e.g., Park et al., 2009) which, at least in the case of Kīlauea, are spatially coincident with interpreted highly porous (40%) magma reservoirs that erupt magmas with deformed olivine phenocrysts (Pietruszka et al., 2015(Pietruszka et al., , 2018;;Wieser et al., 2020).These high velocity zones are typically surrounded by low velocity zones which tend to be collocated with crustal seismicity.Leahy et al. (2010) analyzed P-to-S receiver functions (RFs) on both onshore and offshore seismometers and interpreted other discontinuities within the crust that may act as barriers to rising melt, such as the interface between the aerial and subaerial volcanic loads, the base of the oceanic crust, and the base of magmatic underplating, which seems to be present below all of the Hawaiian Islands (Leahy et al., 2010).Despite this complexity, most studies generally agree that the Moho lies between 12 and 18 km BSL and roughly correlates with topography (Doran & Laske, 2021;Hill, 1969;Leahy et al., 2010;Li et al., 1992;Zucca et al., 1982).
Below the Moho, and particularly in the uppermost lithospheric mantle, it has proven to be somewhat challenging to track the magma plumbing system below Hawaiʻi.Low shear-wave velocities between 20 and 40 km BSL are observed east, northeast, and southwest of Hawaiʻi, with the largest anomalies implying ∼2%-3% melt in the east (Doran & Laske, 2021).Relatively deep seismicity is observed ∼30-50 km depth near the southwestern coast of Hawaiʻi which has been interpreted as large-scale tectonic faulting in the mantle lithosphere primarily due to stresses imposed by the overlying volcanic load with possible stress contributions from magmatism (Wolfe et al., 2003).This faulting has also been hypothesized to promote lateral transportation of magma through the lithosphere (Wech & Thelen, 2015).In addition to this deep seismicity, high rates of tremor have been identified between 36 and 43 km in depth in a relatively small area slightly west of Kīlauea's southwest rift zone (Figure 1b [pink box], Wilding et al., 2023).This location is interpreted as a sill complex potentially representing Mauna Loa and Kīlauea's shared magma supply from the upper mantle (Wilding et al., 2023) and perhaps indicating where the deeper lithospheric mantle plumbing system bifurcates.This interpretation, however, is tenuous given geochemical signatures seem to require separation of their sources below the lithosphere (Abouchami et al., 2005;Frey & Rhodes, 1992;Weis et al., 2011).At greater depths, the lithosphere-asthenosphere boundary (LAB), as constrained by S-to-P RFs, varies from about 75 km beneath the islands to ∼93 km beneath the surrounding oceanic basin (Rychert et al., 2013).The LAB is underlain by a positive amplitude anomaly that deepens from ∼125 km beneath the islands to ∼155 km about 100 km west of Hawaiʻi, interpreted to result from a +100°C temperature anomaly associated with the position of the hotspot plume.This agrees with the location of low velocities in the upper mantle based on teleseismic P-wave and earthquake surface wave tomography (Figure 1b, Laske et al., 2011;Wolfe et al., 2011), although the exact position of these slow velocities vary somewhat depending on methodology (e.g., Wolfe et al., 2009).These low velocities extend down into the mantle transition zone (MTZ), where a relatively thin MTZ is observed in the area (Cao et al., 2011).The thin MTZ and 10.1029/2023JB027029 3 of 13 low seismic velocities are interpreted to result from the penetration, or at least the transfer of heat, of the mantle plume through the MTZ from the lower mantle.The significant lateral migration of hotspot material from the lower mantle to the surface is attributed to rising melts exploiting pre-existing lithospheric weaknesses, likely due to fracturing resulting from lithospheric flexure due to the volcanic load (Hieronymus & Bercovici, 1999;Rychert et al., 2013).
Clearly, the linkage between deep mantle "plume" and the Hawaiian hotspot is geometrically complex.In this study, we aim to constrain the upper lithospheric discontinuity structure below Hawaiʻi down to 50 km depth using P-to-S wave RFs.Ultimately, we aim to improve our understanding of the seismic structure of the region, and interpret our seismic images in the context of seismicity and geochemistry to provide insight into the lithospheric magma plumbing system beneath Hawaiʻi.

Data and Receiver Functions
The RF technique allows for the recovery of Earth's discontinuity structure below a seismic station by isolating converted phases on a seismogram through deconvolution of different components of ground motion (Langston, 1979).In this study, we investigate the crustal and uppermost mantle structure beneath the island of Hawaiʻi using P-wave radial RFs at 47 stations.We used Mw > 5.5 events between 30 and 95 epicentral degrees 2,071) while the green dots/bars indicate the proportion of those events remaining (total: 1,130) after controlling for data quality.Distribution is dominated by events from the west and southwest, but all backazimuths have a similar proportion of accepted to rejected events.(b) Map of seismic stations used in this study are colored by network with major volcanic edifices labeled and their rifts outlined in white and associated faulting outlined in black (U.S. Geological Survey, Quaternary fault and fold database for the United States, accessed 17 April 2023, at: https://www.usgs.gov/natural-hazards/earthquake-hazards/faults).The pink box in the south represents the bounds of the Pāhala sill complex located between ∼36 and 43 km BSL (Wilding et al., 2023).(Inset) A regional mapview of the Hawaiian Islands with locations of low P-wave velocities (<−0.25% from the depth slice average) for depths of 100, 300, and 400 km (Wolfe et al., 2011).10.1029/2023JB027029 4 of 13 from our stations so that waveform complexities due to triplication associated with the MTZ velocity structure and core diffraction can be avoided.We filtered the raw data from 0.04 to 4 Hz and rotated the data to the ZRT reference frame.For quality control, we first visually inspected the vertical components of each event to ensure a clear P-wave was present at each station.Then, we calculated ∼1.4 and ∼2.5 Hz (2.8 and 5.0 Gaussian alpha parameter) P-wave radial RFs to highlight both lithospheric and intracrustal discontinuities beneath the study area.We used time-domain iterative deconvolution (Ligorría & Ammon, 1999) until a misfit of 0.01% or 800 iterations were reached.Finally, RFs were manually inspected to ensure high signal-to-noise ratios and realistic waveforms.We use data downloaded from the IRIS data center at 47 stations and 2,071 events from August 2004 to September 2021, resulting in a total of 33,709 even-station pairs.After inspection, 1,130 events were considered leading to a total of 3,530 event-station pairs (Figure 1).

Adaptive Common Conversion Point (ACCP) Stacking
Common Conversion Point (CCP) Stacking is a method that maps RFs to a 3D volume of velocity discontinuities throughout an area (Dueker & Sheehan, 1997).In this method, RF amplitudes are traced along their theoretical ray paths according to an a priori P-and S-wave velocity model.Model space is then discretized into spatial bins of a given size, and the corresponding amplitudes of the rays that fall within a bin after converting RF time to depth are averaged to represent the discontinuity structure at that location.In this study, we use a modified version of CCP stacking that allows lateral bin size to adapt to ray density allowing the model to resolve finer detail in regions that are densely sampled by the data, while smoothing the model in areas with low data density (adaptive CCP Stacking, or Adaptive Common Conversion Point [ACCP] Stacking; Delph et al., 2015).We raytrace and convert our radial RFs to depth using a linearly interpolated 1D velocity model derived from a local tomography study (Lin et al., 2014) extended to 100 km (Figure S1 in Supporting Information S1).Bin spacing in our model is 0.05° with a bin height of 1 km and minimum bin width of 0.05°.If <5 rays are present in a given bin, then bins are allowed to expand in increments of 0.05° until five rays are present in a given bin, up to 0.25° width.The amplitude of RFs that fall into the same bin are then stacked to create a 3D model of conversion amplitudes throughout the region.

The Effects of Elevation Corrections in ACCP Stacking
Generally, CCP stacks are created assuming that elevation changes below an array are negligible compared to the vertical dimension of the bins.However, the elevation of Hawaiʻi varies from sea level to ∼4 km over a relatively small distance.Given that our bins have a vertical dimension of 1 km, not accounting for changes in elevation will potentially alter the bin that a particular time segment of a ray will fall in leading to the stacking of data that isn't appropriately mapped to depth.To correct for elevation, the ACCP volume and 1D velocity model of Lin et al. (2014) was extended above sea level to the elevation of individual stations.Ray paths were then traced to depth from the station elevation rounded to the nearest divisible number of bin height (e.g., a station elevation of 0.65 km was rounded to 1 km).
We demonstrate the effect of this elevation correction implementation with a synthetic example containing elevation variations similar to those at Hawaiʻi.We model the vertical and radial response of an incoming teleseismic plane wave to a discontinuity below 21 stations with 10 km spacing for a simple layer-over-halfspace model using the teleseismic plane wave modeling code of Frederiksen and Bostock (2000).We place stations at a range of elevations consistent with the stations in Hawaiʻi using a Gaussian function to represent volcanic island topography and crustal thicknesses that range from 20 km below the lowest elevations to ∼32.8 km below the highest elevations at ∼2.8 km.We account for the geometry (i.e., dip) of the interface in the forward modeling, and use a crustal density, P-wave, and S-wave velocity of 2,800 kg/m 3 , 6.8 km/s, and 3.9 km/s and mantle density, P-wave, and S-wave velocity of 3,300 kg/m 3 , 8.1 km/s, and 4.5 km/s, respectively.Events were created with a ray parameter of 0.06 s/km and a backazimuthal distribution from 0 to 360 with 10° separation.We then create P-wave radial RFs and ACCP stacking using the same approach as described for the real data (time-domain iterative deconvolution, bin widths of 0.05-0.25°and heights of 1 km, and a minimum hit count of five rays per bin).For simplicity, we use a 1D velocity model consisting of a P-wave velocity of 6.8 km/s and S-wave velocity of 3.9 km/s to map the RFs to depth (same as used for the forward modeling) as this will appropriately map conversion from the Moho interface to depth.As seen in Figure 2, not accounting for changes in elevation can cause up to a 4 km error in the imaged depth to the Moho for this model, much higher than the theoretical vertical resolution given the frequency content of our RFs.This elevation correction leads to a more accurate mapping of RFs to depth.

Receiver Functions
In general, RFs throughout the island show a relatively strong conversion at ∼1.5-2 s, consistent with previous studies (Leahy et al., 2010) that estimate a Moho discontinuity around 2 s across the island (e.g., Figure 3; Figures S2-S4 in Supporting Information S1).However, some stations showed particularly complex structure as well, especially near the Kīlauea caldera and ERZ.As should be expected, the seismic structure of volcanic ocean islands can be very complex, resulting from pre-existing boundaries, magma reservoirs, and a complex structure that may include structural (i.e., dipping interfaces or faulting) and/or intrinsic (i.e., mineral-controlled) anisotropy.These effects on RF amplitude and arrivals are covered extensively in multiple studies (Frothingham et al., 2023;Schulte-Pelkum & Mahan, 2014;Schulte-Pelkum et al., 2020).This complexity can be observed at station HV.WRM, directly northwest of Kīlauea, which shows a somewhat consistent arrival around 1.5 s along with significant backazimuthal variations with rays crossing beneath Kīlauea showing large amplitude negative arrivals.ACCP stacking can somewhat mitigate this through "stacking out" complexity if a sufficiently diverse set of rays (i.e., from different backazimuths and/or stations) fall within a bin.Thus, while amplitudes in our individual RFs may be affected by lateral heterogeneities, these will likely result in second order structures compared to the common arrivals observed between stations at depth.

ACCP Stack and Cross Sections
We compute ACCP stacks for 5 Gaussian (5G) (Figure 4 and Figure S6 in Supporting Information S1) and 2.8 Gaussian (2.8G) (Figures S7 and S8 in Supporting Information S1) filtered RFs.Using the 5G ACCP stack, we utilized the open-source software OpendTect (dGB Earth Sciences (2021).https://dgbes.com/software/opendtect) to perform horizon tracking on the highest positive amplitude conversion in our ACCP stack as a proxy for Moho depth (Figure 4).Throughout much of the island, variations in the depth to the continuous, high amplitude horizon range from about 12 to 18 km BSL.However, significant variations in the depth to this conversion are also seen below the active volcanic edifices of Mauna Loa and Kīlauea, likely due to the complex structure of the Hawaiian crust (Figures 4 and 5).The cross sections shown transect, or nearly transect, the major volcanic centers of Hawaiʻi and some of their rifts: Kīlauea, Mauna Loa, Mauna Kea, and Hualālai (Figure 4).
We image a strong, positive arrival beneath Kīlauea ∼6 km BSL, about 5 km shallower than the surrounding positive amplitude horizon, surrounded by weak negative conversions and seismicity (Figures 5b and 5d).This negative conversion and consistent seismicity continues beneath Kīlauea's East Rift, slightly dipping seaward WRM is relatively close, slightly northwest, of Kīlauea while HSSD is distant from any nearby edifice but still close to Mauna Loa's northeast flank.Amplitudes for HSSD are scaled to four times that of WRM. and overlays a relatively strong positive conversion (Figure 5c).Seismicity extends vertically below the mentioned strong, shallow positive conversion near Kīlauea to ∼35 km BSL where it disperses.This deeper seismicity is collocated with weaker, positive amplitude conversions; however, these signals arrive at a similar time as predicted multiples from the high amplitude shallow arrival (Figure S9 in Supporting Information S1), making it difficult to know whether this deeper arrival represents structure.At ∼30 km BSL, we observe a negative conversion to the northwest of this deeper seismicity and the caldera (Figure 5b: ∼125 km along section, D: ∼100 km along section).Beneath the Southwestern Rift of Kīlauea between 35 and 40 km BSL, we find a moderate, negative conversion that spatially correlates with the Pāhala seismic cluster (Figures 5a and 5c) (Wilding et al., 2023).
Mauna Loa has a moderate, positive conversion ∼15 km BSL that disappears southwest of the edifice (not shown in cross section) which is overlaid by two stronger, positive conversion at ∼8 and ∼4 km BSL (Figures 5a and 5d).Seismicity is present along this entire vertical structure but dissipates quickly below the ∼15 km conversion.We observe a weak, negative conversion on Mauna Loa's northwestern rift ∼6 km BSL with little seismicity.Two negative conversions also appear at ∼30 and ∼40 km BSL approximately beneath the edifice and are lacking in seismicity as well.
Mauna Kea is characterized by a strong positive conversion ∼18 km BSL which is accompanied by a small negative conversion at slightly shallower depths (∼14 km BSL) to the southeast which may also be associated with Mauna Loa's northwestern rift.A weak negative conversion ∼30 km BSL underlays the strong positive conversion.We find another negative conversion slightly northwest of Mauna Kea ∼25 km BSL which is underlaid by a stronger negative conversion ∼40 km BSL (Figure 5b).These negative conversions are nearly or completely absent in seismicity.

Positive, High Amplitude Horizon
We expect the general crustal structure of Hawaiʻi to consist of a mafic volcanic load underlain by oceanic crust with possible mafic underplating (Leahy et al., 2010;Li et al., 1992;Ten Brink & Brocher, 1987;Watts et al., 1987;Wölbern et al., 2006).Similarities in the seismic properties between these layers make it difficult to image these inferred structures using our RFs (Fametani et al., 1996) although some interpretations have been made (see Leahy et al. (2010)).Throughout much of Hawaiʻi, we interpret the high amplitude horizon from our ACCP stack as representing the Moho boundary, as it agrees with previous refraction studies (Hill, 1969;Zucca et al., 1982) and other Ps RF studies (Leahy et al., 2010;Li et al., 1992).Excluding Kīlauea and a flank of Mauna Loa, the positive amplitude horizon varies by only a few kilometers.Beneath Kīlauea and Mauna Loa, however, the apparent lack of a Moho conversion could be high velocities within the crust resulting from fractional crystallization of mafic melt creating an ultramafic cumulate material (predominantly olivine) within the magma reservoir.These ultramafic cumulates may possibly extend down to the crust-mantle interface increasing crustal seismic velocities.
In many other seismic studies (local/teleseismic tomography and refraction), high P-wave velocities are observed at relatively shallow depths beneath Kīlauea and Mauna Loa, which are spatially collocated with the high, positive conversions beneath these edifices (∼6 km BSL) in our study (Ellsworth & Koyanagi, 1977;Haslinger et al., 2001;Hill, 1969;Hill & Zucca, 1987;Lin et al., 2014;Okubo et al., 1997;Park et al., 2009;Syracuse et al., 2010;Thurber, 1987;Zucca et al., 1982).The location of these shallow high seismic velocities and high amplitude conversions beneath volcanic edifices also roughly agree with a residual gravity study done by Kauahikaua et al. (2000) indicating some degree of correlation between high density materials and high seismic velocities.The origin of these high velocities has been attributed to extensive differentiation resulting in cumulate formation within a magma reservoir (Ellsworth & Koyanagi, 1977;Hill, 1969;Hill & Zucca, 1987;Lin et al., 2014;Okubo et al., 1997;Park et al., 2009).While the positive amplitudes and the high seismic velocities found beneath Mauna Loa in local tomography are collocated, the strength of the conversion is weaker than at Kīlauea.Low ray density and binning parameters around Mauna Loa (Figures S5, S11, and S14 in Supporting Information S1) in our study may explain this discrepancy.The high velocity zones beneath Kīlauea and Mauna Loa suggest the presence of very shallow, fast velocities that approach the velocity we would expect for ultramafic material similar to the .P-wave tomographic images show very fast velocities (Vp up to 7.5 km/s) at depths of ∼6 km collocated with very shallow strong velocity contrasts in the ACCP images and a vertical trend of seismicity (projected within 0.25° of the cross section) below the Kīlauea.These high velocity zones are surrounded by relatively slow P-wave velocities (Vp ∼6.5 km/s) at similar depths, where we observe velocity decreases in the ACCP images.Pyroxenite + liquid thermobarometry estimates, translated to depth using a P-curve from past seismic and density studies for a given P-T estimate from jadeite contents of the pyroxenes (cyan boxes, Hill & Zucca, 1987;Putirka, 1997), also seem to mostly lie deeper than the very shallow high velocity body beneath Kīlauea.The dark blue box indicates the area of an inferred olivine mush pile based on deformed olivine crysts suggested by Wieser et al. (2020) whose depth is constrained through melt inclusion in olivine crysts translated to depth through a pressure density curve (Tuohy et al., 2016) (See Figure S10 in Supporting Information S1 for other cross sections).
10.1029/2023JB027029 9 of 13 mantle based on (a) their continuity with the high amplitude horizons consistent with the Moho at greater depths elsewhere below the island and (b) the lack of a high amplitude below these shallow discontinuities.
While seismic observations clearly indicate the existence of high velocity material below our volcanic edifices, it remains unclear how the compositional structure and other parameters of these reservoirs would produce the observed seismic structure.Ideally, seismic velocities would be complementary to the geodetic observations that identify two primary magma reservoirs beneath Kīlauea as well as connections to the two rifts from the lower reservoir (Poland et al., 2014).Clinopyroxene thermobarometry indicates that some magma stagnation and evolution occurs near interpreted seismic boundaries (Putirka, 1997; Figure 6).The composition of these reservoirs are constrained mostly by geochemical analysis of whole-rock samples from erupted lavas as well as olivine phenocrysts (Clague et al., 1991;Wieser et al., 2020).Differences in magnesium content between whole rock samples and primary melts indicate a significant volume of olivine in the reservoir.Olivine melts and spinel inclusions also support formation at crustal P-T conditions, so erupted olivine crystals did not originate from the mantle but rather from magma reservoirs through differentiation.Olivine deformation provides further insight to olivine crystal interactions and how these crystals are packed (Wieser et al., 2020) for at least the upper portion of the magma reservoirs.This will have direct implications on seismic velocities (Delph et al., 2021;Paulatto et al., 2022;Takei, 2002;Yoshino et al., 2005).Textural features of erupted olivine phenocrysts, however, are inconsistent with low melt percentages, rather supporting a porosity closer to 40% (Wieser et al., 2020).This would appear as an extreme negative arrival in RF imaging which is not observed, indicating that the region where these olivine phenocrysts are deforming is transient and/or thinner than the theoretical vertical resolution of our RFs (∼1 km).Receiver functions are sensitive to the integrated velocity structure over a given depth range defined by their frequency content, which may be dominated by cumulates rather than thin films of high percentage melt, so this could be consistent with pockets of high percentage of melt if distributed over small vertical distances.The existence of this crystalline ultramafic material would result in velocities much faster than typical crystalline crustal material (Christensen, 1996), even if melt is present (Delph et al., 2021;Paulatto et al., 2022), and could explain the large positive amplitude discontinuity at shallow depths and the observed >7.5 km/s P-wave velocities from local tomography.Indeed, seismic velocities would allow for an integrated basaltic melt percentage of up to 5% if hosted in an ultramafic magmatic reservoir, as it would still be faster than the surrounding crust.These velocity-melt relationships do not take into account anelastic conditions produced by higher temperatures that would be expected here, thus the 5% number represents a likely maximum melt percentage at the resolution of our data (∼1 km).

Relationship Between Structure and Seismicity
Seismicity occurs in three major areas around the island of Hawaiʻi: (a) in regions related to volcanism as part of the Mauna Loa and Kīlauea shallow magmatic systems (Klein et al., 1987), (b) along the proposed decollement that is thought to lie at the interface of the volcanic load and pre-existing oceanic seafloor sediments (Got & Okubo, 2003;Got et al., 1994), and (c) well into the upper mantle on broad tectonic-like faults at ∼40 km BSL with or without contributions from associated magmatic pathways and a dense cluster centered at ∼40 km BSL.(Klein, 2016;Wech & Thelen, 2015;Wilding et al., 2023;Wolfe et al., 2003; Figure 7).Volcanic seismicity lies above, below, and to a lesser extent, around the 6 km high amplitude arrival beneath Kīlauea, down to ∼35 km BSL, while the arrival itself shows relatively little seismic activity.This quiet zone has been inferred to be a ductile, cumulate magmatic reservoir (Clague & Denlinger, 1994;Park et al., 2007; Figure 7).Additional volcanic seismicity is also seen beneath Mauna Loa, though it is unclear if a portion of this seismicity could be attributed to the decollement.
Decollements are likely to exist between the volcanic load and oceanic plate, but the difference in seismic velocity is likely not significant enough, relative to the velocity difference at the crust-mantle boundary, to be obvious in RF imaging.It is, therefore, difficult for us to link the seismicity to the volcanic load and oceanic crust interface.However, by assuming an oceanic crustal thickness of 6-7 km, the decollement would exist 6-7 km above the imaged Moho (Lin & Okubo, 2020; Figure 7).A cluster of earthquakes are observed at approximately that depth, agreeing with previous studies on decollements (Leahy et al., 2010;Wilding et al., 2023), but the cross section running along the ERZ (Figure 5c) is the only cross section that shows a rough tectonic-like spatial trend in seismicity that can be associated with the decollement.
Deep seismicity (∼35 km) is explored by Wilding et al. (2023) who clearly depicts the interconnectedness of deeper seismicity with shallow seismicity.In our ACCP volume, we see a slight negative amplitude conversion around the southwestern rift of Kīlauea at about 40 km BSL (Figures 5a and 5c).This negative conversion is collocated with a zone of high seismicity in the upper mantle (the "Pāhala sill complex"; Wilding et al., 2023) and, thus, may represent melt within the sill complex.Unfortunately, in our study, we were unable to image the connectivity of this deeper sill complex with crustal structure using RFs due to either contamination of multiples from earlier arrivals or the nature of magma transportation between the Pāhala sill complex and the lower reservoir beneath Kīlauea being primarily diking.Teleseismic RFs, approaching the stations at a near vertical angle, are not sensitive to vertical structures that we would expect with diking.Other approaches may be necessary to image these magma pathways in which the ray paths sufficiently cover a wide variety of angles between 15 and 40 km BSL.There also remains the mystery of the Pāhala sill complex's relationship with even deeper structure in the mantle lithosphere and its connection to the supposed location of the mantle plume.

Conclusions
The hotspot volcanic island, Hawaiʻi, has been shown to be complex in its structure, generation, and evolution.We use RFs to highlight the discontinuity structure across the main Hawaiian Island.Elevation-corrected ray Green color represents the lithospheric mantle, dark gray represents underplating, brown is oceanic crust, and lighter gray is the volcanic load.Magma reservoir locations are calculated from geodetic data during pre-and post-eruption time frames (Poland et al., 2014) and shaded red regions are possible pathways for magma propagating through the lithospheric mantle to the crustal.The dark blue box represents deformed olivine crystals (Wieser et al., 2020) with depth constrained by melt inclusions in olivine crysts translated to depth through a pressure density curve (Tuohy et al., 2016).Cyan boxes signify P-T estimates of clinopyroxenes thermobarometers converted to depth using P-curves from past seismic and gravity studies (Hill & Zucca, 1987;Putirka, 1997).The lower question mark represents a possible connection of the Mauna Loa magma supply to the PSC as interpreted by Wilding et al. (2023).Question marks in the crust represent possible olivine cumulate piles beneath other edifices for which our study may not have the resolution to image.paths are used to construct a 3-D discontinuity model through ACCP stacking of P-to-S RFs from which a consistent, sharp discontinuity is observed between 12 and 18 km, roughly following elevation, and is interpreted to represent the crust-mantle transition.A high seismic velocity discontinuity arrival deviates from this trend beneath the Kīlauea edifice and, to a lesser extent, beneath Mauna Loa.These high velocity zones agree with previous seismic studies and likely represent an ultramafic cumulate pile with a low melt percentage (<5%) over ∼1 km.This high velocity zone is surrounded by lower seismic velocities and high rates of seismicity while somewhat lacking in seismicity in the arrival itself.This seismicity then extends vertically down to ∼35 km at which point it trends subhorizontal to the Pāhala sill complex where a small negative amplitude is observed in our ACCP stack.Our RFs successfully capture the discontinuity structure of Hawaiʻi improving our understanding of the upper lithospheric magmatic structure associated with oceanic hotspot volcanoes.The facilities of IRIS Data Services, and specifically the IRIS Data Management Center, were used for access to waveforms, related metadata, and/or derived products used in this study.IRIS Data Services are funded through the Seismological Facilities for the Advancement of Geoscience (SAGE) Award of the National Science Foundation under Cooperative Support Agreement EAR-1851048.Brandon Herr was funded through start-up funds to Jonathan R. Delph through Purdue University.I would like to thank Guoqing Lin and others for their local tomography model and John D. Wilding and others for their relocated earthquake catalog.I would also like to thank the reviewers for their insightful and constructive comments.

Figure 1 .
Figure 1.(a) Rose diagram and map view showing backazimuthal event distribution for this data set.Red dots/bars represent all seismic events considered (total: 2,071) while the green dots/bars indicate the proportion of those events remaining (total: 1,130) after controlling for data quality.Distribution is dominated by events from the west and southwest, but all backazimuths have a similar proportion of accepted to rejected events.(b) Map of seismic stations used in this study are colored by network with major volcanic edifices labeled and their rifts outlined in white and associated faulting outlined in black (U.S. Geological Survey, Quaternary fault and fold database for the United States, accessed 17 April 2023, at: https://www.usgs.gov/natural-hazards/earthquake-hazards/faults).The pink box in the south represents the bounds of the Pāhala sill complex located between ∼36 and 43 km BSL(Wilding et al., 2023).(Inset) A regional mapview of the Hawaiian Islands with locations of low P-wave velocities (<−0.25% from the depth slice average) for depths of 100, 300, and 400 km(Wolfe et al., 2011).

Figure 2 .
Figure 2. (a) Non-elevation corrected Adaptive Common Conversion Point (ACCP) stacking for synthetic teleseismic events recorded on an array of seismic stations (black triangles at 0 km depth).Red amplitudes show the mapped location of the Moho while the black line is the true Moho depth.(b) Same as (a) with the elevation correction incorporated.(c) Calculated differences in the highest amplitude on ACCP stacks for the non-elevation (blue) corrected and elevation corrected (red) ACCP stacking.Zero represents the "true" Moho depth in the model used to create the synthetic data (black dotted line).

Figure 3 .
Figure 3. Stacked receiver functions from stations HSSD and WRM corrected for each of their respective event-station ray parameters and sorted by backazimuth.WRM is relatively close, slightly northwest, of Kīlauea while HSSD is distant from any nearby edifice but still close to Mauna Loa's northeast flank.Amplitudes for HSSD are scaled to four times that of WRM.

Figure 4 .
Figure 4. Mapview of shaded topography, station coverage, the major edifices, and structures of Hawaiʻi.White delineated features represent rift flanks and faults are outlined in black.Cross section lines correspond to Adaptive Common Conversion Point stacks in Figures5 and 6.The colormap represents a minimum curvature map of the horizon tracking below sea level in Opendtect.This surface can also be seen in the purple line in each of the four cross sections (Figure5).

Figure 5 .
Figure 5. Adaptive Common Conversion Point stacks of five Gaussian receiver functions (cross section locations shown in Figure 4).Black dots represent the location of detected earthquakes within 0.25° from the cross section from Wilding et al. (2023), who considered continuous waveforms recorded on the HVO network between 01 November 2018 and 01 May 2022.The purple line in each of the four cross sections represents the tracked high amplitude horizon (from Figure 4).Black triangles at the surface represent seismic stations used in this study.

Figure 6 .
Figure 6.Comparison of local P-wave tomography(Lin et al., 2014) (top)  and our Adaptive Common Conversion Point (ACCP) stacking images (bottom).P-wave tomographic images show very fast velocities (Vp up to 7.5 km/s) at depths of ∼6 km collocated with very shallow strong velocity contrasts in the ACCP images and a vertical trend of seismicity (projected within 0.25° of the cross section) below the Kīlauea.These high velocity zones are surrounded by relatively slow P-wave velocities (Vp ∼6.5 km/s) at similar depths, where we observe velocity decreases in the ACCP images.Pyroxenite + liquid thermobarometry estimates, translated to depth using a P-curve from past seismic and density studies for a given P-T estimate from jadeite contents of the pyroxenes (cyan boxes,Hill & Zucca, 1987;Putirka, 1997), also seem to mostly lie deeper than the very shallow high velocity body beneath Kīlauea.The dark blue box indicates the area of an inferred olivine mush pile based on deformed olivine crysts suggested byWieser et al. (2020) whose depth is constrained through melt inclusion in olivine crysts translated to depth through a pressure density curve(Tuohy et al., 2016) (See FigureS10in Supporting Information S1 for other cross sections).

Figure 7 .
Figure7.Our D to D′ cross section underlaid by a conceptual model of Hawaiʻi corresponding to the same cross section with tomography contours and seismicity (black dots).Seismicity from A to A′ (blue dots) is projected onto the conceptual model to highlight the location of the PSC.Green color represents the lithospheric mantle, dark gray represents underplating, brown is oceanic crust, and lighter gray is the volcanic load.Magma reservoir locations are calculated from geodetic data during pre-and post-eruption time frames(Poland et al., 2014) and shaded red regions are possible pathways for magma propagating through the lithospheric mantle to the crustal.The dark blue box represents deformed olivine crystals(Wieser et al., 2020) with depth constrained by melt inclusions in olivine crysts translated to depth through a pressure density curve(Tuohy et al., 2016).Cyan boxes signify P-T estimates of clinopyroxenes thermobarometers converted to depth using P-curves from past seismic and gravity studies(Hill & Zucca, 1987;Putirka, 1997).The lower question mark represents a possible connection of the Mauna Loa magma supply to the PSC as interpreted byWilding et al. (2023).Question marks in the crust represent possible olivine cumulate piles beneath other edifices for which our study may not have the resolution to image.