The arc-scale spatial distribution of volcano erosion implies coupled magmatism and regional climate in the Cascades arc, United States

O’Hara, D.; Karlstrom, L.

doi:10.3389/feart.2023.1150760

ORIGINAL RESEARCH article

Front. Earth Sci., 19 June 2023

Sec. Volcanology

Volume 11 - 2023 | https://doi.org/10.3389/feart.2023.1150760

The arc-scale spatial distribution of volcano erosion implies coupled magmatism and regional climate in the Cascades arc, United States

D. O’Hara¹*

L. Karlstrom²

¹Department of Geography, Vrije Universiteit Brussel, Brussels, Belgium
²Department of Earth Sciences, University of Oregon, Eugene, OR, United States

The morphology and distribution of volcanic edifices in volcanic terrains encodes the structure and evolution of underlying magma transport as well as surface processes that shape landforms. How magmatic construction and erosion interact on long timescales to sculpt these landscapes, however, remains poorly resolved. In the Cascades arc, distributed volcanic edifices mirror long-wavelength topography associated with underlying crustal magmatism and define the regional drainage divide. The resulting strong along- and across-arc modern precipitation gradients and extensive glaciation provide a natural laboratory for climate-volcano interactions. Here, we use 1,658 volcanic edifice boundaries to quantify volcano morphology at the arc-scale, and reconstruct primary edifice volumes to create first-order estimations of Cascades erosion throughout the Quaternary. Across-arc asymmetry in eroded volumes, mirroring similarly asymmetric spatial distribution of volcanism, suggests a coupling between magmatism and climate in which construction of topography enhances erosion by orographic precipitation and glaciers on million-year timescales. We demonstrate with a coupled landscape evolution and crustal stress model that mountain building associated with magmatism and subsequent orographically-induced erosion can redistribute surface loads and direct subsequent time-averaged magma ascent. This two-way coupling can thus contribute to Myr-scale spatial migration of volcanism observed in the Cascades and other arcs globally.

1 Introduction

Volcanic topography represents the interface between crustal magma transport and surface environment. Within volcanic arcs, the distributions of volcanoes and their morphological differences characterize the spatial and temporal evolution of both volcanism and erosion (Völker et al., 2011; Karátson et al., 2012). The spatial distribution, composition, and extruded volumes (e.g., Sherrod and Smith, 1990; Pearce and Peate, 1995; Martin et al., 2004; Till et al., 2019) of volcanic vents and associated edifices reflect the structure of the underlying magma plumbing system (e.g., Vogt, 1974; Hieronymus and Bercovici, 2001), while edifice topography records time-variable volcanic landscape construction and erosion (Karátson et al., 2012; Grosse et al., 2020).

Volcanoes present an attractive target for estimating erosion due to their relative consistency as conical constructional landforms. As volcanoes are pervasive on Earth and other planets (Head and Wilson, 2022), these landforms provide a template for isolating the impacts of climate or primary composition on landscape form. Over the scale of a volcanic arc, rates of edifice erosion provide insight into the topographic evolution and weathering of volcanic terrain. Previous work used large, eroded edifices as passive indicators of past climate conditions (e.g., Karátson et al., 2012) or explored the effects of volcanic arc construction and erosion on global climate (e.g., Lee et al., 2015). However, with the exception of feedbacks on shorter timescales (e.g., glacial unloading; Jull and McKenzie, 1996), possible feedbacks between volcanic topography, magma transport, and spatial climate gradients have not been thoroughly explored. Furthermore, although global compilations of edifice geometries exist (Grosse et al., 2014; Grosse and Kervyn, 2018; Karlstrom et al., 2018), few studies have attempted regional completeness in landform catalogs (e.g., Bablon et al., 2020; Grosse et al., 2020; Bertin et al., 2023).

In this study, we analyze the spatial distribution of edifice morphologies from volcanoes that have erupted in the U.S. portion of the Cascades arc over the past 2.6 Ma using 10 m-resolution Digital Elevation Models (DEMs) and volcano topographic boundaries (O’Hara et al., 2020). We then use two simplified reconstruction methods to estimate lower and upper bounds for edifice eroded volumes, which approximate erosion rates along the arc front. Observed erosion patterns imply coupling between Quaternary Cascades erosion, climate, and volcanic edifice construction. Based on these data we propose a long-term volcano-climate feedback in which mountain building associated with magma transport and eruptions drives orographic asymmetry in erosion, which redistributes topographic surface loads that guide volcanic vent locations. The observed relationships between where volcanoes occur and how they erode thus implies a mechanism for enhancing volcanic arc front migration away from the trench in time. We demonstrate the mechanical feasibility of this feedback using a magmatic landscape evolution model under Cascades-like conditions, where deviatoric stresses in the subsurface from evolving topography guide volcanism.

2 The Cascades arc

The U.S. portion of the Cascades arc spans ∼1,000 km from northern California to Washington (40°N–49°N), and is associated with the subduction of the Juan de Fuca plate under the North American plate (Hildreth, 2007). Quaternary volcanism throughout the arc has culminated in ∼3,000 known vents (Ramsey and Siebert, 2017; O’Hara et al., 2020), spatially arranged in clusters around major volcanoes (defining ‘volcanic centers’; Hildreth, 2007; Karlstrom et al., 2015), or in distributed vent fields throughout the arc (e.g., Bacon et al., 1997; Fleck et al., 2014; Germa et al., 2019). These vents are associated with four dominant edifice classifications (Figures 1A,B)–scoria cones (number of vents n = 2030), domes (n = 357), shield volcanoes (n = 242), and composite volcanoes (n = 206, including large ‘stratovolcanoes’ representing the highest regional topography).

FIGURE 1

FIGURE 1. (A,B): Cascades arc topography and volcanic edifices. Black dots represent all edifice locations, black-dashed line represents the regional drainage divide, and colored symbols distinguish age of last eruption for (A) scoria cones and domes, and (B) shield and composite volcanoes. Regional divide derived by filtering topography over 10 km wavelength and combining drainage basin divides. Abbreviated epoch ages are H.: Holocene, P.: Pleistocene (general), L.P.: Late Pleistocene, M.P.: Middle Pleistocene, E.P.: Early Pleistocene. (C) Topography filtered over a 24 km wavelength (approximately two times the mean stratovolcano diameter) to remove edifice topography. (D) Edifice volume-weighted Gaussian kernel densities ( $λ$ ; see Supplemental text) using a 12 km bandwidth (Connor et al., 2019) to demonstrate edifice distributions. (E) Along-divide distributions of topography (black colors) and edifice volume-weighted Gaussian kernel densities from (D) (red colors) in 10 km swaths across the divide, a proxy for the spatial distribution and volume of volcanic output. Solid lines represent mean values. Black-dashed line represents filtered topography from (C).

Cascades arc topography creates a roughly north-south trending regional drainage divide (Figures 1A–D, black-dashed lines) whose position varies east-west along its length by ∼100 km. As shown through volume-weighted Gaussian kernel densities of edifices (Supplemental text; Connor et al., 2019) within the Cascades (Figure 1D), this divide largely collocates with the arc-front, where the majority of volcanism is focused (Figure 1E). Longer wavelength topography (filtered at 24 km to remove volcanic edifices in Figure 1C) mirrors edifice distributions, as observed at other arcs globally (George et al., 2016). Clockwise tectonic block rotation present in the study area since at least the early Miocene implies extension in the south and compression in the north of the arc, as well as an apparent longitudinal migration of the arc both north and south of the rotation pole (Wells et al., 1998; Wells and McCaffrey, 2013). This kinematic framework is consistent with along-strike segmentation of Quaternary volcanic vents (Guffanti and Weaver, 1988; Wells et al., 1998). Topography mirrors diverse crustal geophysical signals associated with current magmatism in the Cascades arc (e.g., isostatic residual gravity, seismic tomography, heat flux, crustal rotation; Till et al., 2019; O’Hara et al., 2020), so mountain building likely reflects crustal thickening associated with magmatic intrusions that define long wavelength topography (Perkins et al., 2016), as well as background tectonic strain field.

The spatial trend of long wavelength topography and volcano distributions coincides with a significant east-west precipitation gradient (Fick and Hijmans, 2017), with drier conditions and lower erosion rates towards the east (Reiners et al., 2003), as well as a north-south gradient in temperatures and precipitation. Furthermore, repeated cycles of glaciation throughout the Quaternary have been a major component of erosion within the Cascades, especially in the Washington Cascades (Mitchell and Montgomery, 2006). This glacial influence prominently includes Cascades volcanoes, with glaciation being documented along the Cascades arc front as far south as ∼42°N latitude (e.g., Batchelor et al., 2019) and many stratovolcanoes still exhibiting alpine glaciation (Post et al., 1971).

3 Methods

Ramsey and Siebert (2017) compiled several decades of previous field reports to catalog ∼3,000 Quaternary vents in the U.S. Cascades, including the locations, morphologic classification of the associated edifice, bulk composition, and epoch age of last eruption of each vent. O’Hara et al. (2020) utilized this database to delineate 2,105 edifice boundaries from topography using a semi-autonomous algorithm.

3.1 Edifice morphology

Using 10 m-resolution DEMs from the National Elevation Dataset (U.S. Geological Survey, 2013), we isolate edifices by clipping topography around the boundaries previously defined by O’Hara et al. (2020). Afterwards, we collect common metrics associated with bulk edifice morphology (Grosse et al., 2014) for 1,658 edifices (Figure 2A; Supplementary Table S1). These include the planform area ( $A$ ) of the edifice boundary and the overall relief of the edifice ( $h$ ; elevation difference between the lowest boundary and the peak). We also measure three slope metrics of the edifice — the mean local slope (mean slope of each pixel within the boundary), the mean landform slope (mean slope generated from the edifice peak to each boundary point), and the conical slope calculated as $h / \sqrt (A / π)$ . The other 447 edifices either had collapsed calderas (e.g., Mt. Mazama, Newberry) or did not have a clearly-defined, single boundary (e.g., clusters of monogenetic cones that were defined by one boundary), and were thus excluded from the analysis.

FIGURE 2

FIGURE 2. Description of morphology metrics collected for each edifice, using Mt. Adams (WA) for example. (A) Edifice relief, basal planform area, and slope of the entire landform; colors correspond to local topographic slope values. (B) Eroded volume estimate methodologies. Red text and inset describe contour-based eroded volume estimate using a convex-hull. Blue text and shading describe the ray-based eroded volume estimate.

3.2 Eroded volumes

Previous studies have reconstructed edifices and estimated eroded volumes by fitting general volcano shapes over intact, non-eroded surfaces (i.e., planèzes) (e.g., Mitchell and Lofi, 2008; Lahitte et al., 2012; Karátson et al., 2016; Dibacto et al., 2020). However, such methods require extensive knowledge of the edifice’s volcanic history, which is not known for every edifice in the Cascades. Furthermore, our analysis incorporates edifices with sizes that span a range of magnitudes over an entire arc, as opposed to volcano-specific scales explored previously. We thus estimate eroded volumes of the edifices using two simplified methods that attempt to correct for edifice incision and flank mass removal (Figure 2B).

The first method reconstructs the non-eroded edifice shape using a convex-hull approach. In this method, the edifice is divided into a series of elevation contours. Along each contour, the footprint of the edifice is expanded by smoothing concave sections of a polygon while keeping the outermost points stationary, resulting in a ‘convex hull’. Afterwards, the reconstructed edifice is generated by interpolating the new polygons for each contour as a surface. This method relies on closed elevation contours. As edifice topographic boundaries are often irregular and span multiple elevations (e.g., Bohnenstiehl et al., 2012), this method expands higher-elevation portions of the boundary vertically-downward to encapsulate lower edifice flanks with a closed contour.

The second method reconstructs the edifice as an approximate cone-like structure. First, the edifice’s volumetric centroid is determined in map view (Lerner et al., 2020). A series of rays are then projected from the centroid to the edifice boundary in 10° azimuthal bins. Along each ray, the location of highest topography is determined and a straight line is generated between the highest topography and boundary. This line is assumed to be the pre-eroded edifice flank. These flank points are then interpolated to derive the reconstructed edifice.

From both reconstruction methods, the eroded volume ( $V_{e r o d e}$ ) is calculated by subtracting the current topography from the reconstructed edifice and integrating over positive values of the resulting grid. Calculating $V_{e r o d e}$ in this manner, as opposed to subtracting the volumes of the current and reconstructed DEMs, allows us to spatially differentiate and ignore sections that do not fit within the framework of our algorithms (e.g., parasitic cones on the flanks of larger edifices). For each edifice, we also calculate eroded volume as a percentage ( $V_{%}$ ) of the initial landform as $V_{%} = V_{e r o d e} / (V_{t o p o} + V_{e r o d e}) * 100$ , where $V_{t o p o}$ is the volume of current topography. Error analysis using a landscape evolution model for edifice erosion suggests that the ray-based approach is most accurate to determine eroded volumes, but that accuracy decreases with edifice degradation (Supplementary Figure S2). Assuming the edifices within our database are less than 50% degraded, our analysis suggests a maximum error of ∼10% for the ray-based reconstruction and ∼30% for the convex-hull method.

Of the 1,658 edifices analyzed for morphology, we estimate eroded volumes of 1,646 edifices using our automated reconstruction algorithms. The other 12 edifices that could not be adequately processed for eroded volumes had complex topography or planform areas less than 0.08 km² (Figure 4A) such that grid resolution likely impacted the ability of our algorithms to reconstruct topography.

The methodology presented here and in O’Hara et al. (2020) provides a semi-automated workflow for quantifying erosion on volcanic edifices. This approach is not without drawbacks, as edifice identification algorithms often generate conservative boundaries that exclude lower edifice flanks (van Wees et al., 2021; Bertin et al., 2023). As discussed in the Supplemental text, our edifice reconstruction methods assume a single edifice landform that spans variably protracted construction phases and that erosion dominantly occurs post-edifice construction. These methods cannot reconstruct the evolution of total edifice relief, and do not distinguish post-eruptive erosion (e.g., glacial carving, channel incision) from co-eruptive mass wasting (removal of flank material during an eruption; e.g., 1980 Mt. St. Helens eruption). Figure 2 and Supplementary Figure S1 presents an example of the reconstruction methods and eroded volumes for Mt. Adams (WA). These algorithms generate edifice reconstructions that infill erosive features (i.e., valleys) from the current topography.

4 Results

4.1 Edifice morphology

Scoria cones and domes have the lowest mean planform area (∼0.5 km² and 1.0 km², respectively) and relief (∼10 m and 180 m, respectively), as expected for monogenetic landforms (Karlstrom et al., 2018) (Table 1; Figure 3). In comparison, shield and composite volcanoes have mean planform areas of ∼10 km² and 9 km², respectively, and relief of ∼390 m and 330 m, respectively. Of all of the edifice classifications, composite volcanoes have largest range of geometries, reflecting the variability of this classification to span from small composite cones to the well-known, major stratovolcanoes (e.g., Mt. St. Helens, Mt. Shasta, Mt. Rainier; Figures 3D,H).

TABLE 1

TABLE 1. Edifice morphology metrics and eroded volume estimates, arranged by classification. A is planform area, h is edifice relief, and S signifies slope. Edifice erosion rates are calculated by dividing eroded volumes by edifice planform area and mean age for each epoch designation (Holocene/Late Pleistocene = 0.05 Myr; Middle Pleistocene = 0.95 Myr; Early Pleistocene = 2.2 Myr); rates are then averaged over each morphologic classification.

FIGURE 3

FIGURE 3. Edifice geometry Probability Density Functions (PDF) separated by morphologic type. (A–D): Edifice planform area. (E–H): Edifice relief. Solid lines represent mean geometry values. Gray bars are this study; red colors are previously-reported geometries for scoria cones (Karlstrom et al., 2018), shield volcanoes (Grosse and Kervyn, 2018), and composite volcanoes (Grosse et al., 2014). Blue lines in (D, H) represent edifice area and relief ranges of Cascades stratovolcanoes specifically.

Although mean planform areas differ ∼1–2 orders of magnitude between monogenetic and polygenetic edifices, relief is similar. This leads to differences in the mean slopes for each edifice classification (Table 1). Both landform slope and conical slope values (Supplementary Figure S3) follow expectations for each classification — shield volcanoes have the lowest mean landform slope (∼17°), domes have the highest mean landform slope (∼25°), and scoria cones and composite volcanoes have moderate mean landform slopes (∼18° and 21°, respectively). Differences likely reflect average eruption style and composition of each landform class (noting that intrusions also contribute to edifice form; Cosburn and Roy, 2020). Mean local slopes follow a similar trend, but are 3°–7° lower than the mean landform slopes, likely reflecting shallower gradients near edifice summits (e.g., craters) and lower flanks (Grosse et al., 2014; McGuire et al., 2014).

4.2 Edifice erosion

We find that between 50 $\pm$ 15–245 $\pm$ 24 km³ of material has eroded from Cascades edifices in Quaternary time (not including known eruption-related mass removal; e.g., Mt. St. Helens, Mt. Mazama) (Table 1; Figures 4A–C), with the convex-hull reconstruction algorithm associated with the lower bound and ray-based reconstruction associated with the upper bound. Within each edifice morphologic classification, erosion rates spans 5–8 orders of magnitude with significant variability between classifications (Supplementary Table S1). As suggested by Karlstrom et al. (2018), edifice planform areas are well-correlated with eroded volumes through power-law relationships (Figure 4A). Given that edifice areas range 3–5 orders of magnitude within each classification (Figures 3A–D), this large range of eroded edifice volumes is thus not surprising, and demonstrates the impact of landform size on eroded material. Across classifications, eroded volumes increase from scoria cones and domes to shield and composite volcanoes (Figure 4B). Furthermore, edifices associated with low and intermediate silica content (basalt, andesite) have higher eroded volumes than more evolved compositions (dacite, rhyolite) (Figure 4C).

FIGURE 4

FIGURE 4. Estimated edifice eroded volumes. (A) Relationship between edifice area and eroded volume. Colors correspond to reconstruction method. Colored-dashed lines are power-law relationships. Black dots represent edifice areas that could not be adequately reconstructed. (B,C): Estimated eroded volume ranges (bars) and means (horizontal lines) for each (B) morphologic type and (C) composition. (D) Bar graph showing the cumulative convex-hull (blue bars) and ray-based (red bars) eroded volume estimates as a percentage of reconstructed edifice volume for each morphologic type. Bars represent the total volume budget (sum of all eroded volumes divided by the sum of all reconstructed volumes across the classifications) and lines represent mean erosion percentages (average values of eroded volume divided by reconstructed volume for each edifice).

The apparent increase in eroded volumes between edifice morphologies likely reflects primary edifice volumes rather than erosional efficiency. Figure 4D displays eroded volumes as percentages of reconstructed edifice volumes for each classification to remove possible bias associated with landform volume. Total eroded volume percentages range 2%–4% and 7%–15% for the convex-hull and ray-based reconstruction methods, respectively, with composite volcanoes having the highest total percentage. Mean eroded volume percentages vary between 0.8%–1.2% and 3.6%–7.1%. For both reconstruction methods, domes experience the lowest mean erosion percentages and shield volcanoes have the highest.

4.3 Spatial distribution of edifice erosion

A regional Quaternary landform and edifice erosion catalog allows us to evaluate the spatial distribution of volcanic output and erosion within the Cascades at 0.5° latitudinal bins. Figures 5A,B show the along-arc distributions of estimated eroded volumes and total reconstructed edifice volumes. Figure 5C shows eroded volume as percentages of reconstructed edifice volume.

FIGURE 5

FIGURE 5. Along-arc spatial distribution of edifice eroded and total reconstructed volumes, in 0.5° latitudinal bins. (A) Convex-hull (black bars) and ray-based (red bars) estimated eroded volumes. (B) Total erupted volume based on edifice topography. (C) Edifice eroded volume as a percentage of the total erupted volume. (D) Edifice eroded volume percentages across the regional drainage divide. Black triangles (convex-hull estimations) and red stars (ray-based estimations) represent distribution values with the addition of erosion from eruptive process for Mt. Mazama and Mt. St. Helens, with black-dashed arrow ambiguity related to Mt. Mazama eroded volume (Bacon and Lanphere, 2006) (see Supplemental text). Blue shading in (D) represents along-arc precipitation rate ranges across the divide, blue lines are mean values.

The distribution of volume-normalized edifice erosion (Figure 5C) suggests apparent segmentation along the arc. The most northerly section (48–49°N) encompasses Mt. Baker and Glacier Peak, and has eroded volume percentages ranging 2%–7%. The segment between 46–47.5°N spans from Mt. Rainier southward to the Columbia River Gorge, and contains the highest erosion percentages, up to ∼22%. South of the Columbia River Gorge, erosion percentages exhibit a general southward decrease in values, from ∼19% at 45°N to ∼2% at 40°N. Although erosion percentages vary between bins, a noticeable break in this trend occurs at 42.5°N, the latitude corresponding to the location of Mt. Mazama (i.e., Crater Lake) which has a collapsed caldera and was thus excluded from both the morphology and erosion analysis. The break in erosion percentage at this location thus likely reflects the lack of eroded volume estimates, rather than additional segmentation (Guffanti and Weaver, 1988; Schmidt et al., 2008).

Figure 5D presents the along-arc distributions of eroded volume percentages split across the topographic divide. The same along-arc trends in eroded volume percentages can be observed, but there is a distinct east-west asymmetry, with eroded volume percentages 1.1x–7.9x higher on the western part of the arc than on the eastern.

4.4 Post-eruption edifice erosion rates

Although well-dated eruptive histories exist for only a fraction of Cascades volcanoes (Hildreth, 2007; Schmidt and Grunder, 2009; Wall et al., 2019), the last eruption of each vent in our Quaternary vent database is known at the epoch scale. Combining these data with our estimated eroded volumes allows for a first-order estimation of edifice erosion rates post-dating eruptive activity over the last ∼2.6 Myr (Supplemental text). Complex and perhaps non-monotonic construction-erosion histories of polygenetic edifices means that erosion occurring between active episodes cannot be accounted for. For example, flank instabilities in the form of lahars and dome collapses (Miller, 1980) help shape the morphology of edifices and derived basal contours, altering derived eroded volumes and erosion rate. However, with some exception (Grosse et al., 2022), individual eruptions or episodic flank collapses do not typically fully resurface edifices (Bacon, 2008). We therefore derive time-averaged edifice erosion rates by dividing eroded volume by the edifice’s planform area and mean age of the epoch associated with the most recent eruption (2.2 Myr for the Early Pleistocene, 0.95 Myr for the Middle Pleistocene, and 0.05 Myr for the Late Pleistocene and Holocene). Estimating erosion rates in this manner is an oversimplification, but appropriate given the highly-variable age constraints and the regional, arc-wide scope. Furthermore, we note that, incomplete geologic data notwithstanding, our erosion rates are comparable to more detailed studies in other continental arc settings (10^–3–10^–1 mm/yr; Karátson et al., 2012; Dibacto et al., 2020).

Supplementary Figures S4A, B shows the range and mean erosion rates for each edifice classification and vent composition. Across all volcanoes, mean erosion rates range in magnitude between 10^–3–10¹ mm/yr. Similar to the eroded volume estimates (Figure 4B), there is an overall increase in mean erosion rates from monogenetic (10^–3–10^–2 mm/yr) to polygenetic edifices (10^–2–10^–1 mm/yr), topped by composite volcanoes. Separating erosion rates by composition suggest basaltic and andesitic edifices have higher mean rates (10^–3–10^–1 mm/yr) compared to dacitic and rhyolitic edifices (10^–3–10^–2 mm/yr). Being normalized by edifice basal area, the increase in erosion rate between classifications likely reflects thresholds in erosion processes such as drainage area (e.g., Montgomery and Dietrich, 1988); the differences in erosion rate between compositions may reflect variable primary rock reactivity, fracture density, or topographic roughness.

Calculating eroded volumes and time-averaged erosion rates based on the last known eruption of a volcano is most accurate for monogenetic edifices, which often form either through a single (sometimes years-long) event, or through multiple events on 0.1–1 kyr timescales (e.g., Inbar et al., 1994; Báez et al., 2017; Smith and Németh, 2017). Because monogenetic edifices are numerous throughout the arc, this provides an opportunity to quantify erosion through time at slightly higher temporal resolution. First, considering average monogenetic eroded volumes associated with each epoch (Supplementary Figure S4C), we find that eroded volumes follow the expected trend of increasing with time. However, this trend is non-monotonic, with a larger increase in eroded volume between the Holocene/Late Pleistocene and Middle Pleistocene than between the Middle and Late Pleistocene. Mean time-averaged monogenetic edifice erosion rates from each epoch (Supplementary Figure S4D) mirror this trend and show decreasing rates through time, with younger landforms (Holocene and Late Pleistocene) having time-averaged rates of 10^–2–10^–1 mm/yr and older edifices (Late Pleistocene) having time-averaged rates of 10^–4–10^–2 mm/yr. These trends match previous studies that suggest volcano erosion follow logarithmic relationships (e.g., McGuire et al., 2014; Bablon et al., 2020).

5 Discussion

5.1 Global- and regional-scale volcano landform distributions

A discrepancy between edifice morphologies occurs when comparing distributions of edifice planform areas and relief from our study to previously-published global datasets (Figure 3). Scoria cone and shield volcano relief distributions are similar to other studies (Grosse and Kervyn, 2018; Karlstrom et al., 2018); however, overall edifice planform area and relief distributions are one to two orders of magnitude smaller than those previously reported for other locations. Two possibilities exist to explain these differences between datasets. The first is that the edifice boundary identification algorithm used by O’Hara et al. (2020) (based on the Modified Basal Outlining Algorithm; Bohnenstiehl et al., 2012) creates more conservative boundaries than those used in the previous studies (e.g., NETVOLC; Euillades et al., 2013), which would thus lower both edifice planform areas and relief. However, a simple comparison between major stratovolcano geometries between this study and those reported by Grosse et al. (2014) shows that although our boundaries create lower morphometry values, they are generally within the same magnitude (Supplementary Table S2).

Rather, we argue that the differences between datasets highlight an important overestimation of global volcanic landform catalogs. Previously-reported global scoria cone morphologies compiled older studies from 7 volcanic fields that spanned a range of settings (Karlstrom et al., 2018), while shield and composite volcano morphology compilations (Grosse et al., 2014; Grosse and Kervyn, 2018) were dominantly derived from the Smithsonian Global Volcanism Program (Global Volcanism Program, 2013). This overestimation is most apparent with shield and composite volcano planform areas, where only the largest edifices of our database overlap with previously-reported data. This suggests these global studies are incomplete at the arc-scale, giving a biased interpretation towards larger edifices.

5.2 Regional climate coupled with volcanic topography

In the Pacific Northwest U.S., a latitudinal gradient in modern (30-year average) precipitation (Fick and Hijmans, 2017), correlated with spatial variations in average winter temperature (Ersek et al., 2012), aligns with the direction of the Cascades (Figure 5D; 6C). Precipitation also varies strongly in an across-arc direction due to the orographic effect, with paleoenvironment and exhumation analyses suggesting these patterns likely existed over 10-Myr timescales (Reiners et al., 2003; Bershaw et al., 2019). These precipitation gradients are reflected in reconstructed edifice volumes and erosional percentages. For example, Figures 6A,B plots convex-hull reconstructed volumes and eroded volume percentages, respectively, relative to precipitation rate as binned averages for each classification (Supplementary Figure S6 provides same figure for ray-based eroded volumes), demonstrating that reconstructed edifice volume and erosion correlate with precipitation (particularly for polygenetic edifices). Relationships between topography, erosion, volcanic vent distributions, and precipitation also occur spatially (Figure 7).

FIGURE 6

FIGURE 6. Comparison of edifice and climate data. (A,B): Convex-hull based (A) reconstructed volumes and (B) eroded volume percentages as function of mean precipitation rate for each edifice morphologic type, given as averaged bins in mean precipitation intervals of 0.25 m/yr. (C) Relationships between polygenetic volcano convex-hull eroded volume percentages (circles) and mean precipitation rate (background grid) in geographical reference frame with horizontal exaggeration. Black line represents regional drainage divide, light-gray shading and red lines denote estimated glacial extent during the Last Glacial Maximum (Ehlers et al., 2011; Watt et al., 2013; Crandell, 2015; Batchelor et al., 2019) (Supplementary Figure S8).

FIGURE 7

FIGURE 7. (A) Spatial relationships between volcanic edifice, eroded volumes, and climate data, showing polygenetic volcano convex-hull eroded volume percentages (colored circles) and mean precipitation rates (background grid) projected relative to distance from regional drainage divide. Red outlines denote major Cascades stratovolcanoes. (B) Comparison between precipitation, maximum elevation, and volume-weighted Gaussian kernel densities in 20 km bins east (red triangles) and west (blue squares) of the divide. Colored regions denote data extent between east and west. (C) Binned across-arc comparison of (lower panel) polygenetic vent frequency, (middle panel) polygenetic edifice convex-hull eroded volume percentages, and (upper panel) mean elevation (black line) and precipitation rate (blue line), relative to the arc axis (black-dashed line).

Our data show that significant asymmetry in edifice erosion across the divide mirrors a ∼2-fold difference in mean precipitation values (Figure 5D, blue lines) within the southern and central sections of the arc, suggesting a long-term influence of orographic precipitation on edifice morphology. Because the topographic divide collocates with the volcanic arc axis in all but the northern Cascades (Figures 1C,D), across-arc relationships between arc topography, precipitation, and erosion can be analyzed spatially by considering values relative to the divide (Figure 7).

Projecting precipitation rates, polygenetic edifice locations, and convex-hull eroded volume percentages in map-view relative to the divide (Figure 7A) demonstrates that the accumulated mass of volcanic edifices spatially controls precipitation patterns. This relationship is further highlighted by plotting maximum elevations, mean precipitation rates, and edifice volume-weighted Gaussian kernel densities in 20 km bins along both the eastern and western portions of the arc (Figure 7B). Our data suggest that topography, volcanic construction, and precipitation are well-correlated on the eastern part of the divide — higher volcanic vent densities correspond with both higher precipitation and topography. However, these relationships are not as clear in the west, with higher vent densities corresponding to higher elevations, but over a range of precipitation values.

Taking across-arc bins of composite volcano vent and eroded volume distributions, Figure 7C better resolves these relationships by showing that although the majority of Cascades composite volcanoes cluster around the divide (bottom panel, black bars) and correlate with topography (upper panel, black line), the highest average eroded volume percentages along the arc front are spatially offset towards the west of the divide by ∼20 km (middle panel, red bars). This offset also corresponds with the highest mean precipitation rates within the arc (upper panel, blue line), consistent with predictions of orographic precipitation (Smith et al., 2005) based on regional average atmospheric circulation patterns and topography (Supplemental text).

Edifice eroded volume percentages also increase with latitude (Figure 5C; Supplementary Figure S5), correlating with both higher precipitation rates as well as glaciation throughout the Quaternary (e.g., Batchelor et al., 2019) (Figure 5D; 6C). As evidenced by glacially-sculpted modern stratovolcano topography (e.g., Mt. Rainier, Mt. Washington, Sisters Volcano Complex), scattered volcanic landforms reflecting direct ice-magma interactions (e.g., tuyas in the central Oregon Cascades; Schmidt and Grunder, 2009) and moraine deposits (O’Connor et al., 2001), a significant fraction of volcano erosion during the Quaternary is glacial in origin. Last Glacial Maximum (LGM; ∼26.5–19 ka; Watt et al., 2013) ice is well constrained to have extended through the Washington Cascades with up to ∼1 km thickness (Siegert, 2001), but farther south is less certain. Figure 6C shows the cumulative extent of glaciation from previous reconstructions of the LGM (Supplemental text), indicating variably-focused glaciation on the arc axis through southern Oregon.

Similar to other arcs (e.g., George et al., 2016; Perkins et al., 2016), volcanic edifices correlate with long-wavelength topography in the Cascades (Figure 1D). This is expected based on magmatic intrusive-to-extrusive ratios (e.g., White et al., 2006) and implied by modern geophysical anomaly patterns (O’Hara et al., 2020). Given the relationships between edifice spatial distributions and crustal magmatic structure, as well as between eroded volume percentages and climate, we view Cascades edifices as proxy for the regional-scale interaction between erosion patterns and magmatism. Correlations shown in Figures 5D, 7 suggest two important relationships between arc-scale climate and volcanism. First, mountain building associated with magmatism controls across-arc precipitation gradients. Second, this across-arc precipitation gradient creates spatial asymmetry in surface erosion between the windward and leeward sides of the arc. Both of these relationships have been well-documented in other environments, including arcs (England and Molnar, 1990a; Ferrier et al., 2013; Muller et al., 2022), and can feed back mechanically on the shallow magma plumbing system, coupling climate and magmatism on long timescales.

Although it is not clear whether glacial erosion is distributed evenly across the arc axis, this suggests another mechanism for volcano-climate coupling, in which volcanism promotes glaciation that persists through glacial cycles by building high topography that can sustain glaciers. A type of glacier ‘buzz-saw’ effect (Mitchell and Montgomery, 2006), erosion by ice could limit the height of volcanic arc topography and account for increasing erosion rates at higher latitudes.

5.3 Long-term evolution of volcanic arc landscapes

Magmatic emplacement, edifice construction, and volcanic deposition proceed on a spectrum of time and spatial scales that differ from mountain building in purely tectonic settings (e.g., England and Molnar, 1990b). Indeed, volcanic terrains exhibit an inverted relationship between regional uplift and erosion (where intrusion-related uplift can be more localized than erosion; e.g., O’Hara et al., 2019), as well as anomalous river network structures (Sweeney and Roering, 2017) and hydrology (Jefferson et al., 2006; 2010), compared with other active landscapes.

Although current geomorphic transport laws cannot fully encapsulate all of the processes that occur within volcanic environments (e.g., Karátson et al., 2010), strong correlations observed in the Cascades between arc topography, precipitation, and edifice erosion suggest climate-driven long-term feedbacks that could help frame arc landscape evolution. Such feedbacks have been suggested in other orographic settings (e.g., Beaumont et al., 2001; Roe et al., 2003; Pelletier et al., 2010; Stalder et al., 2020; Val and Willenbring, 2022), but have not been explored where landscape evolution is largely driven by magmatic and volcanic processes. Paleoelevation data and the record of ancestral eruptions indicate high volcanic topography in the Cascades since 20–30 Ma (Bershaw et al., 2019), such that feedbacks may be generalizable to similarly long timescales.

5.3.1 Climate feedbacks on arc volcanism

In convergent settings, topography has been argued to encode the interplay between climate and tectonics (e.g., Beaumont et al., 2001; Willett, 2010; Fox et al., 2014; Schildgen and Hoke, 2018). In arcs, topographic form depends on the rate of magmatism relative to regional climatic and tectonic factors (Whipple and Meade, 2006; Perkins et al., 2016; Stalder et al., 2020). From the spatial distributions of both volcanic edifices and eroded edifice volumes, we hypothesize that long-term climate-volcanism coupling may influence patterns of magmatism. We focus on arc front migration as one possible outcome of this coupling (Karlstrom et al., 2014; Muller et al., 2022).

Volcanic arcs exhibit large-scale spatial migration of volcanic output on ∼10 Myr timescales, a phenomenon that has been attributed to slab dynamics (e.g., Dickinson and Snyder, 1978; Karátson and Timár, 2005; Phipps Morgan et al., 2008; Dibacto et al., 2020), magmatic crustal thickening (e.g., Karlstrom et al., 2014), and overriding plate tectonics (e.g., Wells and McCaffrey, 2013). If arc migration reflects pathways for magma ascent that are modulated by the stress state of the crust, it follows that time-evolving surface loads could influence crustal stresses and thus also play a role in directing arc migration.

Orographic precipitation driven by arc mountain building focuses erosion on the windward side of the arc front, and provides a means for such long-term modulation of surface loads. Models of landscape evolution in response to regional climate gradients suggest that non-uniform precipitation causes regional-scale divide migration towards drier regions (e.g., Roe et al., 2003; Goren et al., 2014; Muller et al., 2022) in orogenic settings. Similarly, topographic loading has been shown to influence underlying shallow crustal stresses and modify dike trajectories that guide volcanic vent locations and their migration towards the load through time (e.g., Kervyn et al., 2009; Rivalta et al., 2019).

We test the influence of evolving, climate-driven topographic loading on magmatic shallow dike propagation using a simplified coupled landscape evolution–crustal stress model. Below, we describe the model setup, results, and interpretation in regard to Cascades edifice erosion data.

5.3.2 Landscape evolution–crustal stress model

We develop a landscape evolution model to show the phenomenology of divide migration driven by magmatism and orographic precipitation. Topographic evolution or magmatic patterns of the Cascades are not addressed specifically, although we do use average characteristics of Cascades topography and modern precipitation in our calculations. Furthermore, we do not attempt to explicitly differentiate between surface growth driven by intrusive magma uplift and mantling by extrusive volcanic deposits, focusing of magma ascent within the crust (e.g., Karlstrom et al., 2009), or the time-variable nature of erodibility and hydrology in volcanic terrains (e.g., Jefferson et al., 2010). Our kinematic approach shows that this coupling between climate-driven topographic evolution and crustal stresses that guide shallow magma ascent is possible for a system resembling the Cascades.

The model (described fully in the Supplemental text) couples a 1D landscape evolution model (O’Hara et al., 2019) to a 2D crustal stress model (Styron and Hetland, 2015), with the mechanical feasibility of individual model components being demonstrated in prior studies (Maccaferri et al., 2010; Goren et al., 2014; Styron and Hetland, 2015; O’Hara et al., 2019). The Cascades are idealized as paired drainage basins separated by a hillslope domain. Model elevations ( $z$ ) as a function of distance (x) are calculated through time ( $t$ ) as:

\frac{\partial z (x, t)}{\partial t} = U_{t e c t o n i c} (x, t) + U_{m a g m a t i c} (x, t) - E (x, t) (1)

where $U_{t e c t o n i c}$ is a regional uplift rate (assumed tectonic in origin), $U_{m a g m a t i c}$ is a local surface growth rate that parameterizes both intrusion and eruptions (termed here as magmatic uplift), and $E$ is spatially-variable erosion containing contributions from bedrock river incision (the ‘stream power law’; Whipple and Tucker, 1999) and linear hillslope diffusion (Culling, 1960). Time-evolving topography is used to estimate crustal stresses ( $σ$ ) based on Boussinesq solutions for a vertical load distribution imposed on a homogenous elastic halfspace (Styron and Hetland, 2015).

We use a model domain of 80 km (approximate width of the Cascades arc front) and create an initial model topography as the steady-state solution with uniform precipitation and $U_{t e c t o n i c}$ = 0.1 mm/yr (Moon et al., 2015, based on the north Cascades). In all model phases, $U_{t e c t o n i c}$ is constant in time and space. Following O’Hara et al. (2019), the first model phase constructs volcanic topography using a parabolic uplift function with a scale length of 20 km and maximum $U_{m a g m a t i c}$ of 0.7 mm/yr (mean current stratovolcano relief divided by 2.6 Myr). This uplift rate averages faster uplift rates at volcanic centers (Cascades stratovolcanoes have lifetimes of ∼1 Myr or less; Hildreth, 2007) with surrounding terrain where magmatic uplift is likely slower and intrusion-dominated. As the distribution of Cascades polygenetic edifices is skewed towards the east of the divide (Figure 7C), we center the zone of magmatic uplift 7 km east of the initial divide. This phase is run for 2.6 Myr with a time step of 10 kyr. Afterwards, the second phase is generated by keeping $U_{m a g m a t i c}$ active and incorporating a non-uniform precipitation gradient using mean normalized across-arc precipitation values presented in Figure 7C. This phase is run in 100 kyr time steps until steady-state topography is re-established. Model parameters are discussed in the supplemental text, and the complete model is shown in a Supplementary Video.

5.3.3 Model results

The first phase of our model (Figure 8A; Supplementary Video V1) demonstrates how topography is constructed by a generalized zone of uplift associated with magmatism. These simulations assume parameters defining ‘Regime 3′in O’Hara et al. (2019), in which the uplift perturbation is significant compared to background uplift and the available stream power of the fluvial network. This uplift asymmetrically alters topography to generate steeper slopes (and thus higher erosion rates) windward of the perturbation, thus initiating divide migration towards the locus of magmatism (cf. O’Hara et al., 2019). The second phase (Figure 8B) shows how across-arc gradients in precipitation further drive divide migration. As has been demonstrated by previous studies, the divide migrates away from the wetter portion of the landscape (Bonnet, 2009; Goren et al., 2014; Muller et al., 2022). As localized magmatic uplift is still active in this phase, our model demonstrates that precipitation gradient effects can migrate the divide further eastward than the main locus of volcanic construction.

FIGURE 8

FIGURE 8. Landscape evolution model demonstrating coupling between climate, magmatism, and topographic evolution. Phases are (A) magma-driven uplift builds topography and (B) climate gradient causes erosional asymmetry that generates divide migration, modulating subsurface deviatoric stress trajectories. Upper panels show kinematic forcing functions (uplift and precipitation rates). Middle panels show instantaneous (solid black line) and reference (black-dashed line) topography of each stage. Lower panels show crustal stresses associated with topographic loading (Styron and Hetland, 2015). Colors are magnitude of maximum compression ( $σ_{1}$ ), short black lines are principal stress directions, and solid black lines show stress trajectories as proxy for dike pathways. Maximum topographic stress magnitudes in this model are in excess of 10 MPa within the uppermost km. Peak magma-driven uplift in (A) is offset 7 km east of the divide to simulate across-arc polygenetic edifice distributions in the Cascades (Figure 7C, bottom panel).

In both phases, we see that evolving topography modulates stress orientations and magnitudes in the shallow crust. Stress trajectories perpendicular to the least compressive principle stress axis imply dike propagation pathways (e.g., Anderson, 1951; Muller et al., 2001; Maccaferri et al., 2010; Rivalta et al., 2015) and our model shows that magma ascent is focused towards intrusive uplift (Figure 8A). As climate gradients drive asymmetric erosion and alter the magnitude of divide migration, the effect is mirrored in the subsurface (Figure 8B). This kinematic model demonstrates that climate-driven arc migration over 10s of km can occur on 10s of Myr timescale, similar to rates observed globally under other subduction conditions (Karlstrom et al., 2014). This mechanism may occur over even shorter timescales as well (Muller et al., 2022), likely depending on the relative magnitudes (and spatial distribution) of magmatic uplift and asymmetric orographic erosion.

We note that glaciation has also been proposed to influence magmatism independently, modulating surface loads on Milankovitch timescales to influence magma transport, crustal storage, and mantle melting (Jull and McKenzie, 1996; Jellinek et al., 2004; Hooper et al., 2011). In the Cascades, such inferences are complicated by an incomplete volcanic record (Watt et al., 2013) and low rates of magmatism compared to highly magmatic terrains such as Iceland (Jull and McKenzie, 1996) that exemplify glacio-volcanic coupling. However, because Quaternary glacial cycles in the Cascades are focused at high-elevations along the arc front, time-variable ice loading and erosion should further enhance the relationship between climate-driven asymmetric erosion, crustal loading, and magmatism proposed here.

5.4 How does climate influence volcanic arc migration?

Patterns of volcanic migration in the Cascades are largely consistent with tectonic block rotation of the Pacific Northwest relative to stable North America, which generates an extensional regime in the south and compressional regime in the north of the arc (Wells et al., 1998). Although geologically-inferred displacement vectors generally correspond with modern GPS velocities (McCaffrey et al., 2013), differences have been attributed to variable roll-back of the Juan de Fuca slab under North America and decoupling of upper plate deformation from underlying mantle (Wells and McCaffrey, 2013; Humphreys and Grunder, 2022). Superimposed on this and other geodynamical mechanisms would be contributions from climate to influence crustal magmatism and arc migration, but evidence for such surface influence does exist.

For example, eastward regional drainage divide migration, opposite the direction implied by tectonic rotation, has been suggested in the north Cascades through capture of the formerly eastward-flowing Skagit river by the windward side of the range (Simon-Labric et al., 2014). Muller et al. (2022), based on apparent westward vent migration of recent (<1 Ma) eruptive centers in this region, suggest that magma ascent pathways shift windward (west) in response to leeward (east) drainage divide migration implied by thermochronology. On the other hand, arc-scale Quaternary volcano data compiled here suggest the opposite direction of arc-front migration. Vents are concentrated leeward of the drainage divide on average (Figure 7C, black bars), while edifice erosion is consistently larger on the windward side (Figure 7C, red bars). Because volcanic topography largely defines the drainage divide in the U.S. Cascades (Figure 1), our data imply regional leeward migration of this divide.

On the basis of transcrustal visco-elasto-plastic numerical models of magma ascent, Muller et al. (2022) propose that deep and asymmetric faults flanking the arc are activated through asymmetric surface erosion, which can facilitate magma migration windward as well as leeward. Our model neglects deep structural control and instead approaches magma ascent as through isolated dikes, suggesting that leeward migration of topography alters stresses in the shallow crust to promote migration of magmatism in the same direction (Figure 8). The two approaches are different in significant ways.

Migration of vents is inferred based on asymmetry in visco-elasto-plastic fault patterns (assuming incompressible elasticity) in Muller et al. (2022), which arise from deep crustal deformation forced by asymmetric surface loading. This numerical approach is common to large-scale geodynamic models (e.g.; Keller et al., 2013). In contrast, our model imposes magmatic uplift kinematically, and vent migration is inferred on the basis of (compressible) elastic stress trajectories. This numerical approach is similar to geodetic modeling of active volcanoes (e.g., Anderson et al., 2015). Similarly, Muller et al. (2022) incorporated landscape evolution models that are initialized in a non-steady-state configuration and impose orographic precipitation as a function of surface slope alone. Our landscape evolution approach in contrast isolates parameters kinematically from initial steady-state solutions, which facilitates analysis but is challenging to apply to real landscapes. We see great potential in future studies to understand the impact of these approaches on the emergent feedbacks between magmatism and climate.

Quaternary arc-scale data assembled here seem to favor a leeward vent migration model. For example, despite having relief similar to the northern portion of the Cascades, across-arc precipitation gradients are steeper in the central Cascades. Within a 40-km band centered on the arc axis, precipitation decreases leeward across-arc by a factor of ∼2.5 in the central portion of the Cascades, spatially correlating to the highest amount of arc migration since the Early Miocene (Supplementary Figure S9). In comparison, the northern Cascades have ∼1.5x decrease in across-arc precipitation and lowest distance of arc migration. If modern precipitation rates are proxy for longer-term climate signals, this along-arc relationship between precipitation gradient strength and arc migration distance further supports our proposed feedbacks between climate and magmatism. However, further work may better establish a robust empirical basis for widespread windward migration and asymmetric faulting as proposed by Muller et al. (2022), and seen through the apparent migration of individual eruptive centers on ∼100 s kyr timescales (Hildreth, 2007).

6 Conclusion

Volcanic arcs represent natural laboratories to explore the long-term relationships between magmatism, climate, and topography. Using a database of Quaternary volcanism and satellite-derived modern topography in the Cascades, we analyzed edifice geometry distributions and spatial patterns of erosion throughout the 1000-km U.S. portion of the arc. Relating edifice spatial distributions and first-order eroded volumes to modern precipitation rates and past glacial extents, we find that although Quaternary volcanism aggregately builds topography and cluster around the main arc axis, erosion is highest ∼20 km west of the axis and follows a strong across-arc precipitation gradient.

Expanding on previous studies of divide migration and topographic influence of upper crustal magmatism, we propose a new coupling between climate and magmatism. Our results suggest that arc construction by volcanic activity builds a topographic barrier that generates precipitation gradients and drive divide migration by erosional asymmetry, which can feed back into the magma plumbing system through crustal loading. Numerical landscape evolution–crustal stress modeling demonstrates the feasibility of such feedbacks. Our results thus provide new avenues of exploration and testable hypotheses to better understand the long-term evolution of volcanic arcs in relation to climate factors.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material. The code for the landscape evolution—crustal stress model is available on Github (https://github.com/danjohara/Landscape_Evolution_Crustal_Stress). Further inquiries can be directed to the corresponding author.

Author contributions

DO’H conducted the data collection and analyses, generated figures, and wrote the manuscript. LK advised on the project direction, edited the manuscript, and secured funds for the project. All authors contributed to the article and approved the submitted version.

Funding

NSF CAREER Grant 1848554, awarded to LK.

Acknowledgments

The authors wish to acknowledge discussion with J.J. Roering, N. Klema, and manuscript comments from M.A. Jellinek. Reviews by Karoly Nemeth and two reviewers were valuable and helped focus our presentation.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2023.1150760/full#supplementary-material

References

Anderson, K. R., Poland, M. P., Johnson, J. H., and Miklius, A. (2015). Episodic deflation–inflation events at kı lauea volcano and implications for the shallow magma system. Geophys. Monogr. Ser. 208, 229–250. doi:10.1002/9781118872079.ch11