Time-evolving surface and subsurface signatures of Quaternary volcanism in the Cascades arc

Diversity in the spacing, volume, and morphology of arc volcanoes (e.g., Tamura et al., 2002; George et al., 2016) implies diversity in underlying crustal magmatism. Mapping active structures through the crust to connect volcanism with deeper magmatic processes remains an outstanding challenge. Here, we combined a database of mapped Quaternary vents, surface topography, and diverse geophysical data sets within the Cascades arc (western North America) to probe relations between volcanism and underlying crustal structure. Building on prior efforts to synthesize geophysical (e.g., Weaver et al., 1989; Wells et al., 1998; Till et al., 2019) and geologic (e.g., Guffanti and Weaver, 1988; Hildreth, 2007) data in the Cascades arc, we analyzed (1) arc-scale relations among geophysical data sets associated with magmatism; (2) the extent to which volcanoes match geophysical subsurface magmatic signatures; and (3) temporal variations in these relations during the Quaternary.

Volcanism in the north-south–trending Cascades arc is associated with eastward subduction of the Juan de Fuca plate under the North American plate (Fig. 1A). We focused on the United States Cascades (∼40°N–49°N). Quaternary volcanism consists of notable long-lived (∼300–600 k.y.; Calvert, 2019) stratovolcanoes aligned parallel to the trench, as well as voluminous off-axis volcanic fields encompassing thousands of vents extending as far as ∼50–150 km normal to the trench (Guffanti and Weaver, 1988; Hildreth, 2007). Although clockwise rotation of western Oregon has migrated the arc on ∼10 m.y. time scales (Wells et al., 1998; du Bray and John, 2011), previous work has not documented consistent Quaternary vent migration (Hildreth, 2007).

Mapping of the Cascades has revealed Quaternary volcanic products that span the range of common edifice types and compositions observed on Earth (e.g., Sherrod and Smith, 2000; Hildreth et al., 2012). Ramsey and Siebert (2017) compiled a database containing 2999 vent locations (Fig. 1A; see the Supplemental Material¹), along with associated morphological classification and epoch age of the most recent eruption (Holocene, 0–0.01 Ma; late Pleistocene, 0.01–0.1 Ma; middle Pleistocene, 0.1–1.8 Ma; early Pleistocene, 1.8–2.6 Ma).

We compiled the following Cascades geophysical data sets that document crustal attributes at <∼20 km depth and that may constrain magma structure.

(1) Isostatic residual gravity anomaly data provide a depth-integrated measure of upper-crustal rock density (Blakely et al., 1997), correcting observed gravity for topography and compensating the crustal root (Simpson et al., 1986). We did not seek signatures of magmatic crustal thickening (Karlstrom et al., 2014).

(2) Seismic tomography involves a combination of rock composition, temperature, and fluid content (Zhao et al., 1992). Several tomographic models exist for the Cascades; we used 10 s and 15 s period phase velocity anomalies (ΔV_ph) from a surface-wave model based on both onshore and offshore data, which are sensitive to upper-crustal structures (Janiszewski et al., 2019).

(3) Heat-flow measurements reflect conductive and advective heat transport in the upper few kilometers of crust, with lateral heat advection by groundwater over a scale of tens of kilometers (Ingebritsen and Mariner, 2010).

(4) Crustal rotation rates derived by regional GPS velocity field measurements record interseismic surface motions (McCaffrey et al., 2013) and approximate large-scale rotation rates over the past ∼16 m.y. (Wells and McCaffrey, 2013).

To identify surficial signatures of volcanism in the Cascades, we used 10-m-resolution National Elevation Data set digital elevation models (DEMs; U.S. Geological Survey, 2013) to determine topographic extents of volcanic edifices. Generally, edifices are positive topographic structures associated with vents, they are semicircular in plan view, and they have slopes higher than surrounding topography (including any satellite vents). Volcano topographic boundaries and volumes were determined by the modified basal outlining algorithm (MBOA; Bohnenstiehl et al., 2012). We ignored topography associated with dispersed tephra and lava flows, although deposit volume likely scales with edifice volume. Furthermore, we did not account for buried vents, nor did we account for syn- or postconstruction erosion; however, distributions of edifice volumes were similar for all epochs in our data set (Fig. S5 in the Supplemental Material), suggesting erosion does not bias our results.

We augmented the MBOA with a procedure that uses regional slope hypsometry to isolate volcanoes, calculate edifice volumes as the integral of bounded topography, and subtract small structures such as parasitic cones from underlying edifices (see the Supplemental Material, and Fig. S1). Assuming detailed geologic mapping is generally more accurate, we used estimates from Hildreth (2007) and Bacon and Lanphere (2006) for major stratovolcano volumes in the subsequent analysis. Our topographic volumes generally compared well with Hildreth (2007) and the volcano DEM analysis by Grosse et al. (2014), although some significant differences exist for volcanoes with complex shapes (see the Supplemental Material).

We analyzed cinder cones, domes, shield volcanoes, and composite volcanoes, giving a total of 2835 analyzed vents. Of these, we determined boundaries for 2105 vents. The remaining edifices have morphologies that are not easily distinguishable from surrounding topography. In this study, we assumed that these unidentified vents have volumes equal to our average volumes for each morphologic type, as determined by the MBOA.

We calculated a total minimum Quaternary edifice volume of ∼2730 km³, implying a minimum extrusion rate of ∼1.05 km³/km/m.y. for the ∼1000 km length of the study area and 2.6 m.y. of the Quaternary (Fig. S4, Tables S2 and S3). Figures 1B and 1C show the spatial distribution of edifice numbers, volumes, and extrusion rates. Figure 1D shows the cumulative arc-scale volume of edifices by epoch.

Our estimated volumes are nearly identical to the ∼2570 km³ estimated by Sherrod and Smith (1990) for the U.S. Cascades, which included dispersed deposits. Hildreth (2007) updated this, estimating Quaternary erupted volume of the entire Cascades to be ∼6400 km³. If Hildreth’s estimate is correct, then current edifice volumes account for ∼50% of total Cascades output. Although glacial erosion is variably significant (Hildreth, 2007), if we assume that missing volume comes mostly from deposits, total extruded volumes are roughly twice the volume of edifices alone.

We linearly interpolated regional geophysical data to a common 25 × 25 km² resolution grid (Fig. S6), approximately equal to twice the median diameter determined for major Cascades volcanoes (see the Supplemental Material). Sensitivity tests (Fig. S13) showed that grid resolution did not affect our results.

We then performed a series of correlation calculations. We first assessed structures not explicitly associated with vents (Fig. 2A) by considering regional gridded data sets alone, over the area plotted in Figure 1A. Next, we linearly interpolated gridded data to the analyzed 2835 vent locations to identify structures underneath volcanoes (Fig. 2B). Finally, we subdivided vent data into epochs (Figs. 2C and 2D; see Supplemental Material). Because the vent database records only the most recent eruption for an edifice, we limited temporal analysis to monogenetic vents to mitigate bias from long-lived volcanoes. Temporal variation between Holocene and early Pleistocene vents is discussed here to illustrate variations in vent distribution. All epochs are presented in the Supplemental Material; Figures S8–S12 show the variation between data sets as a sequence of biplots along with associated best-fitting trend lines.

Both the number and magnitude of correlations substantially increase among regional data sets when interpolated to vents compared to regional grids alone (Figs. 2A and 2B). The most significant relations are consistent with a magmatic origin (Fig. 2E). For example, magma-driven temperature or melt anomalies should contribute to lower seismic velocities and isostatic residual gravity while increasing surface heat flux. These relations are all observed (Fig. 2B), and additionally corresponded to higher vent density, larger edifices, and increased elevations (e.g., Cao et al., 2016; Deng et al., 2017). Rotation rate is uncorrelated to other regional gridded data, yet it strongly covaries when interpolated to vents, suggesting magmatic influence on crustal deformation near volcanoes. We therefore interpret correlations between vent density and these geophysical data sets as defining magmatic structures in the upper crust linked to volcanic expression.

Correlation magnitudes are generally higher for Holocene versus early Pleistocene vents (Figs. 2C and 2D; Fig. S7). In spite of coarse temporal resolution, this decrease suggests that older edifices no longer overlie active magmatic structures, especially considering that monogenetic vents are most numerous in the earliest epochs (Table S3).

Figure 2 indicates that surface-subsurface correlations exist at volcanic edifices, but it does not reveal arc-scale patterns. We examined this spatial structure using independent metrics of surface and subsurface data. We assessed surface data with a volume-weighted Gaussian kernel function λ(x,y) (see the Supplemental Material) that measures spatial vent density and edifice volumes as a probability density function (e.g., Connor et al., 2019). We also measured the extent to which subsurface data provide a coherent indication of magmatic structure, but models that relate data physically (e.g., gravitational admittance or Nafe-Drake curve) are not similarly comparable. Therefore, we assessed a relative extent of magmatic influence between data sets with linear bivariate relations.

The correlation of geophysical data as measured by G is largest in central Oregon and generally increased to the south, with more subdued peaks associated with the Caribou (California) and Simcoe (Washington State) volcanic fields, Medicine Lake (California), and Mount Mazama and Mount Hood (Oregon) (Fig. 3A). This pattern is mimicked but more focused in λ (Fig. 3B), in part because weighting vents by volume localizes λ around the large edifices. Broad monogenetic vent fields are also prominent, illustrating the significant distributed volcanism in the central and southern Cascades.

Finally, we note that edifice volumes and vent spatial density distributions covary, both peaking around the Mount Shasta/Medicine Lake latitude. The extent of variation relative to this area thus measures distributed versus focused styles of volcanism. Normalized vent number and edifice volume distributions are plotted in Figures 3C and 3D, along with their difference (β). Positive β implies volumes distributed across more edifices, while negative β indicates volume focused around fewer vents. As expected, volcanic fields such as Caribou, Medicine Lake, and Simcoe are distributed, while areas such as Mount Shasta and Glacier Peak (Washington State) are more focused (Fig. 3D). Vent focusing in areas otherwise dominated by distributed volcanism occurs at Mount Mazama and Newberry volcano (Oregon).

To characterize regional-scale spatial variability in surface volcanism, we compared maximum λ and β among epochs (Fig. 4; where the Holocene was included with the late Pleistocene). These temporal bins are larger than the major Cascades edifice total ages (Calvert, 2019) and so constrain transient patterns of volcanic effusion on million-year time scales.

Figure 4 indicates that arc-scale patterns of volcanism style and magnitude varied throughout the Quaternary. Along-arc patterns of volcanic output are consistent between epochs (Fig. 4A). High λ values cluster in ≤ km-scale areas, highlighting long-lived magmatic centers (Guffanti and Weaver, 1988; Hildreth, 2007). Although lower λ values were found in the early Pleistocene at most of these centers, we could not disentangle true flux variations from vent exposure bias. However, a maximum in λ at the latitude of Lassen Peak (California) in the early Pleistocene likely indicates decreased eruptive output through time in that region.

Along- and across-arc β values hint at changes in volcanic style through time (Figs. 4B and 4C). Vent patterns are not uniform throughout the arc, although a general tendency seems to be northward evolution toward more focusing. Particularly intriguing are two locations where volcanism shifted in style across the arc at the same latitude. Mount Shasta has tended toward focused vents, while Medicine Lake in the rear-arc has become more distributed through time. Exactly the opposite temporal progression was observed ∼300 km north at Newberry volcano and Three Sisters. Although precise dates are lacking, both on- and off-arc axis volcanism may occur simultaneously (Germa et al., 2019).

We speculate that focusing of rising magma is a self re-enforcing process throughout the crust. Radial focusing of vents at Mount Mazama over ∼40 k.y. may have been influenced by thermomechanical feedbacks among volcano loading, pressurized magma storage zones, and rising dikes (Karlstrom et al., 2015). Such organizing processes could operate over length scales of tens of kilometers (Pinel and Jaupart, 2000; Karlstrom et al., 2009), where such vent clustering is observed elsewhere. Tectonic extension, increasing in magnitude south and eastward along the Cascades arc (Guffanti and Weaver, 1988; Schmidt et al., 2008), should promote distributed volcanism and counterbalance focusing. Deeper variations in magma influx to the lower crust (Till et al., 2019) may also influence overlying crustal transport.

We demonstrated the efficacy of remotely deriving edifice volumes from DEMs, as well as arc-scale signatures that relate volcanism to the subsurface. Combining edifice volumes with geophysical inferences of shallow crustal structure, we generated a suite of linear predictors for active magmatic transport pathways under volcanoes along with two metrics that elucidate crustal magmatic structures and spatiotemporal variations in magmatism throughout the arc. The temporal resolution and broad array of subsurface constraints compiled here thus provide a baseline for future efforts to map and model crustal magma transport in the Cascades and other volcanic provinces.

We thank C. Connor, A. Germa, A. Grunder, and M. Guffanti for detailed and constructive reviews. This work was supported by National Science Foundation (NSF) grant GRF-1309047 to O’Hara and NSF CAREER grant 1848554 to Karlstrom. O’Hara completed much of the topographic analysis during a U.S. Geological Survey–NSF GRIP fellowship.