A Multidisciplinary Study on the Crustal Nature of Volcanic Conduits and Magma Reservoirs

Volcanic settings vary widely not only in their eruptive style and products, but in the manner magma travels from deep sources to individual eruptive centers. Imaging these pathways, and their associated crustal reservoirs, provides unique and unprecedented views into these environments. Imaging techniques are varied with the strength of the technique often based on data availability. As such, we focus on two methods—gravity and seismic—in two different settings, each with its own unique volcanic environments, crustal structures, and associated data resources. The f i r s t , t h e Hawaiian Islands, are the most geologically studied hot-spot islands in the world, yet the only largescale compilation of marine and land gravity data is more than 45 years old. We present a new chain-wide gravity compilation allowing us to locate current and former volcanic centers, major rift zones, a previously suggested volcano, and show that volcanoes along the chain are composed of a small proportion of intrusive material (<30% by volume). At the second area, the arc-volcanism of southern Washington, we used ambient seismic noise methods to constrain the crustal pathways of deep-sourced melt to the surface. We image two zones of reduced velocity, one of which correlates with a proposed extensive zone of mid-crustal partial melt which likely supplies evolved magmas to the surrounding volcanoes and vents, including Mounts St. Helens and Adams.


Introduction
The main Hawaiian Islands evolve from active volcanoes on the southeastern end-Mauna Loa, Kilauea and Lo'ihi-to the eroded remnant of Ni'ihau Volcano 600 km to the northwest ( Fig. 1.4a). Progressive cooling and crystallization of magma in crustal reservoirs and surrounding r i f t zones produces rocks rich in olivine (cumulates). These cumulate cores define the long-term average zones where magma resided or transited through, p r i o r to surface eruptions or emplacement in shallow intrusions . Encompassing the cumulate cores are larger zones of dense, dike-rich, intrusions-intrusive complexes-comprising the magmatic plumbing system feeding each volcano. Density contrasts between these features and encompassing lava flows result in positive residual gravity anomalies (Strange et al., 1965;Flinders et al., 2010) and often correlate with fast seismic velocities ). Here, we invert a new chain-wide compilation of land and marine gravity data to estimate the volumes, average densities, and olivine percentages of intrusive complexes and cumulate cores underlying all known volcanoes throughout the Hawaiian Islands.

Data Reduction
Our new compilation is composed of 4820 land-based measurements, including historical data (Strange et al., 1965) and data from recent studies of Kaua'i (Flinders et al., 2010) and the island of Hawai'i . Science and Technology ( Fig. 1.1, inset). These data were filtered to eliminate high-frequency noise, due to changes in survey speed and course, and corrected for crossover errors using x2sys, a part of the Generic Mapping Tools (Wessel and Smith, 1991). The standard deviation of the corrected crossings was 2 mGal.
Complete Bouguer anomalies were produced by removing the effects of topography/bathymetry using a two-part prism-based terrain correction (Fig. 1.2;Flinders et al., 2010). Within 500 km of each measurement, the water column was in-filled with submarine crust ( Watts and Cochran, 1974;Flinders et al., 2010).

3D Inversion
For inversion, the marine residual gravity data were down-sampled onto a 500m cell-spaced geographic grid, resulting in 42,326 measurements. The cell-spacing used in down-sampling the data was approximately equal to half the minimum one-km horizontal wavelength typically resolvable by the marine data. All land residual gravity data were used. The compiled data set was subdivided into three overlapping geographic regions; Kaua'i/Ni'ihau, O'ahu through Maui, and the island of Hawai'i. Data for each region were inverted to produce 3D models of subsurface density contrast using GRAV3D (GRAV3D, 2007). The subsurface was discretized into a set of 3D voxels, each with a constant density contrast relative to 2.7 g/cm 3 , with voxels spanning 1 km in the horizontal and varying in the vertical dimension from 500 m at the surface to 1500 m at depth. The top of the model space was bound by topography/bathymetry while the model basement, at 20 km below sea level, encompassed the base of the thickest Hawaiian volcanic crust (13 km) and 7 km of preexisting oceanic crust (Watts et al., 1985). The density distribution was found by solving an optimization problem of minimizing a model objective function while generating synthetic data that fall within the uncertainty of the observed data (e.g. Flinders et al., 2010). Inversions were subject to the constraint that densities be between 2.0 (wet sand) -3.3 g/cm 3 (olivine). The three individual inversion models were then merged into one chain-wide model ( Fig. 1.4). These inversions provided a low misfit to the residual anomalies, ≈2 mGal rms, and inverting the data with a wide range of initial model parameters verified the returned inversion structure.

Intrusive Complexes and Cumulate Cores
Negative residual gravity ( Fig. 1.3) and low crustal densities (<2.7 g/cm 3 , Fig. 1.4d) are associated with several of the previously mapped debrisavalanche deposits and slumps  cumulate core is likely to be underestimated because it is still an active volcano and the cumulate core is presumably still accumulating.

Rift Zones and Unrecognized Volcanoes
Rift-zone related intrusive complexes are observed beneath the major submarine ridges and tend to be linear features that trace back to individual volcanoes (2.80 g/cm 3 isosurface; Fig. 1 Two unique gravity anomalies/intrusive complexes are seen beneath Ka'ena Ridge, a broad eastern anomaly and a more distinct western anomaly, neither of which appears to extend from the Wai'anae rift zone ( Fig. 1.4a). Sinton et al. (2013) showed that most of Ka'ena Ridge is chemically, structurally and chronologically distinct from Wai'anae Volcano, and likely represents a precursor volcano to the island of O'ahu, consistent with speculations of . Although the broad eastern gravity anomaly has no associated bathymetric high, it is likely associated with this precursor volcano (Ka'ena Volcano), with its older structure possibly buried by younger Wai'anae flows. The western gravity anomaly is located 3 km from the center of the 11-km diameter large cone, which is thought to h ave formed subaerially, given the presence of massive 'a'a flows (Coombs et al., 2004).
Geologic samples acquired during a previous ROV dive, appear to be geochemically similar to the rest of the ridge (Coombs et al., 2004;Sinton et al., 2013). However, given the distinction between this western gravity anomaly and the broad eastern anomaly ( Fig. 1 (Francis et al., 1993;Dzurisin et al., 1984). These estimates were based on geologically short-term observations of the active Kilauea Volcano, specifically heat loss over the Kupaianaha lava pond (1986-1992Francis et al., 1993) and uplift of Kilauea summit (1956Dzurisin et al., 1984). Given our conclusion of primarily extrusive growth, the minimization of the chain-wide residual gravity anomaly by a submarine density of 2.7 g/cm 3 , and that the majority of each volcano is built during it's tholeiitic shield stage (≈95%; Clague and Dalrymple, 1987), we conclude that the majority of each volcano (>70%) is likely comprised of submarine extrusive flows formed during the main shield stage. The disparity between previous estimates of extrusive/intrusive ratios and our own may be due to the limited time (yrs) and localization (Kilauea Volcano) of prior observations or a change in the extrusive/intrusive ratio throughout a volcano's growth.
Voluminous extrusive growth also contrasts with the dominantly intrusive nature of continental volcanism, as well as the formation of normal oceanic crust and flood basalt provinces (e.g., Crisp, 1984;White et al., 2006). In continental settings, magma travels greater path lengths through relatively thick and low-density continental crust and thus is more likely to intrude (White et al., 2006). Oceanic flood basalt provinces also require magma to penetrate thick crust (≈30 km or more; Crisp, 1984), additionally having numerous source dikes that span a wide region, leading to a larger intrusive contribution. At mid-ocean ridges, most of the crust is constructed within a vertical accretion zone, proximal to the ridge axis, and the thickness of the intrusive material is controlled by the local temperature structure (Hooft and Detrick, 1993). In contrast, Hawaiian volcanism originates from magma that penetrates a relatively thin crust <20 km (oceanic crust plus the volcano; Watts et al., 1985). Additionally, magma travels from depth to the near surface primarily through a single central vertical conduit  and typically one or two rift zones. While intrusions occur in the central conduit and rift zones, magma retains sufficient mobility such that the major volume of the volcano is formed from lavas erupted well away from these localized sources.      b For the 2.85 g/cm 3 isosurface, volumes rounded to nearest 50 km 3 ; for 3.00-3.10 g/cm 3 isosurfaces, volumes rounded to nearest 10 km 3 . c Reservoir volume as a percentage of the volcano volume is the ratio of the calculated reservoir volume to the volcano volume. Helens and Adams do? Is there a relationship between the melt pathways and the magma compositions at these three volcanoes? To address these questions, we solved for the crustal shear-wave velocity (Vs) structure of the Mount Rainier area using a novel 3D full-waveform tomographic method 17,19 and all available ambient noise data for the region (see supplementary). We account for phenomena rarely considered but increasingly important at this scale, including; the complex 3D spatial sensitivity of wave propagation, scattering of shortperiod waves by topography, and P/S wave-velocity cross-dependence 18,19,20 .
Despite the geologic complexity of the region, detailed correlations are seen between the velocity structure and previously mapped lithologies (Fig. 2.2a).
Beneath Mount Rainier we observe a 7 by 11 km slow VS zone (< 85% VS_Exp. of 3.64 ± .02 kms -1 , for diorite at 5 km and ~175°C) 21 oriented NW/SE, displaced slightly east of the summit ( Fig. 2.2a, △ MR). The zone extends from 1 to 8 km depth, with an average Vs of 3.0 kms -1 , is aseismic, and coincides with the previously observed slow VP central reservoir 8 . These depths agree with the > 7 km magma source depth, derived from melt inclusion volatile concentrations, in andesites from the 2.2 Ka eruption 6 . No second, mid-crustal, reservoir was found beneath the volcano, and although traces of marginally reduced velocity extend to 16 km, they quickly trend back to the regional basement velocity.
A much larger, NS trending, slow VS zone is located ~15 km west of Mount Rainier, coinciding with the WRSZ (Fig. 2.2a). While the contiguous portion terminates ~ 7km north of the volcano, isolated pockets extend to the shallow reservoir imaged by MT 1,2 . Unlike beneath Mount Rainier, the slow VS zone can be seen well past 16 km depth, forming a gradual extension of a larger mid crustal zone associated with the SWCC (Fig. 2.3). Additionally, an eastward prolongation of the zone connects to the volcano's central reservoir (Fig. 2.2b).
We propose that melt ascends from the deeper crust to within the WRSZ, undergoes crystallization and mixing, before moving laterally eastward into the central reservoir. The slow seismic velocities, aeromagnetic low, and increased seismicity are caused by both the ascending melt, and the large Eocene sediment layer into which it intrudes.
Lateral migration of magma from this off axis-western reservoir to the Numerical modeling supports the conclusion that at subduction zones, the ascending hydrous-basaltic melt phase is emplaced as a succession of sills into the lower crust 4 . These sills generate "hot-zones", where the mantle-derived basalt undergoes various degrees of crystallization to produce residual intermediate-to-silicic H2O-rich melts, as well as cause partial melting of preexisting crustal rocks 4 . As these "hot-zones" evolve, gradients develop in the zone's temperature, pressure, H2O content, and melt fraction, producing the wide range of geochemical variation seen in arc magmas 4 . Further inhomogeneity could also be driven by the varying depth to the subducting slab beneath the "hot-zone" (Fig. 2.1 We can also estimate the melt-fraction of the SWCC by assuming its velocity is a deviation from the eastern regional mid/deep-crustal bedrock. Confining the SWCC to a contiguous body, separate from Mount Rainier's western sediment/melt reservoir, yields an estimated volume of > 28,300 km 3 , with an average reduced VS of 8.4% (Fig. 2.3). Each percent of interstitial melt leads to a decrease in VS of ~1.5 -2% 25 , leading to an average melt-fraction of 5.1 ± 1.3%. At this melt fraction, theoretical temperature estimates 4 are only 50 -75°C above the assumed geothermal gradient 17 , implying a minimal thermal contribution of < 0.5% to the reduced velocity. This volume of melt is more than twice the total eruptive volume of the six main Quaternary eruptive centers from Mounts Hood to Rainier, combined 22 .
We propose that the partial-melt SWCC province is the primary reservoir of of the region's Quaternary arc magmatism in the deep crust, and the crustal expression of subduction-driven mantle melange diapirs 30 . Compositional variations in the volcanoes are, in part, controlled by their location relative to the SWCC, and the subsurface pathways of the ascending hydrous-basaltic melt phase. Mount Rainier stands on the SWCC's diminishing northern margin, which coincides with an abrupt ending of the nearly continuous Mount Hood to Rainier segment of the arc. The next Quaternary cluster to the north (south of Glacier Peak) is separated by a 120 km-long-gap ( Fig. 2.1). We hypothesize that the volcanic segments and gaps may be controlled by SWCC-like crustal "hot-zones" throughout the entire 1100 km-long Cascades arc. Differences in volcanism between the Mount Hood to Rainier segment and the neighboring Garibaldi Volcanic belt to the north and the Oregon Cascades Segment to the south could be partially attributed to the size and emplacement of these crustal "hot zones".
Further imaging of the middle to deep crust across the length of the arc is needed to confirm these crustal provinces. Correlating their locations/extents to Quaternary volcanic centers/intrusions will provide valuable insight on the timescales involved for active volcanism and material cycling in subduction zones.

Methods
Source data is derived from the cross-correlation of ambient noise waveforms, a well proven technique to reconstruct Empirical Green's Function between two seismic-station pairs, i.e. one stations far-field response to an impulsive source at the paired station. While application of ambient noise records to volcanic seismic tomography has become more popular, wave modeling in general has been limited to classical approaches (non-3D/ray-based methods). Ray-based methods assume that the propagation time of a simulated wave from source-to-receiver is only sensitive to velocity perturbations along a line-path connecting the two locations. This assumption does not account for the 3D spatial sensitivity of wave propagation, which is necessary to accurately invert for a complex 3D velocity model. Despite being relatively novel, 3D fullwaveform methods have been shown to yield improved data fits over classical non-3D/ray-based approaches 31,32 .
Our full-wave tomographic method accounts for complex wave propagation using an iterative approach. Each iteration is comprised of (1) finite-difference wave propagation simulations, (2) measurement of phase delays between observed (EGFs) and synthetic waveforms, (3) calculation of 3D P/S-wave sensitivity kernels using the strain Green tensor scattering-integral approach 12 , (4) inversion for velocity perturbations and (5) updating the velocity model. We include a brief discussion on preparation of the EGFs and simulation/inversion parameters, and full details of the method can be found in Gao and Shen (2014).

Empirical Green's Functions
Empirical Green's functions (EGFs) were extracted from short-period

Full-Waveform 3D Simulation and Inversion
Empirical Green's functions were modeled synthetically by simulating the propagation of a Rayleigh surface wave between two station-station pairs using a cartesian version of the 3D nonstaggered-grid finite-difference method of

Delete Earthquakes
Segmentation F.T. Normalization      Shorter-periods are more sensitive to shallow structure, while long-periods have maximum sensitivity more deep within the model.   Sensitivity to small-scale structure (~10 km), is visible to 10 km, but has faded by 25

DATA PROCESSING
km. At these depths we can still recover longer wavelength structure, such as the SWCC. Rainier magma reservoir with and without a neighboring WRSZ reservoir. Smearing is seen between the two reservoirs, indicating we do not have the resolution to comment on their connection from tomography alone.