The Chief Joseph dike swarm of the Columbia River flood basalts, and the legacy data set of William H. Taubeneck Open Access

Dikes represent the primary pathway between magma reservoirs at depth and the upper crust or eruptive centers at the surface. They are the most common mode of magma transport in nearly all volcanic provinces, from oceanic islands (e.g., Krumbholz et al., 2014) to arcs (e.g., Petford et al., 1993) and rifts (e.g., Passarelli et al., 2014). The largest dike swarms exposed on continents are commonly associated with large igneous provinces (LIP). In this style of volcanism, the connection between dikes at depth and lava flows is often poorly exposed (Mège and Korme, 2004; Svensen et al., 2009). As a result, crustal magma transport pathways in these largest terrestrial volcanic systems remain poorly understood. Hence, it is difficult to connect theoretical models of magma flow and dike-related stress changes to observations. Such models are required to assess both the deep mantle origins (e.g., White and McKenzie, 1995) and the surface environmental impacts (e.g., Self, 2006) of these large events.

This study focuses on the dikes that fed the CRBG, the youngest flood basalt province on Earth. The CRBG covers an area of ∼210,000 km² with an estimated volume of ∼210,000 km³ (Fig. 1; Reidel et al., 2013). An estimated 99% of the area was erupted between ca. 17.2 and 15.9 Ma (Jarboe et al., 2008; Barry et al., 2013; Reidel et al., 2013; Camp et al., 2017a; Kasbohm and Schoene, 2018; Cahoon et al., 2020). The Grande Ronde Basalt (GRB), 71% of the CRBG by volume, erupted in ≤400,000 years (Kasbohm and Schoene, 2018). Because most published research has focused exclusively on the extrusive component of the CRBG (e.g., Swanson et al., 1979; Reidel et al., 2013; Camp et al., 2017a, and references therein), connecting surface flows to upper-crustal storage and transport remains an outstanding challenge.

Petrologic mass balance provides one possible approach for estimating intrusive volumes, and because the GRB lavas are notable from other units of the CRBG for their evolved compositions, ranging up to 57% SiO₂, we use this formation for illustration (Hooper, 2000; Reidel and Tolan, 2013). The GRBs have been variously modeled as melts of foundering or recycled mafic crust (Hoshi and Takahashi, 1999; Camp and Hanan, 2008) or as residual liquids produced by differentiation at crustal depths from parents represented by Imnaha Basalt (Wolff et al., 2008; Wolff and Ramos, 2013). In the latter case, a significant addition of mass to the crust in the form of a volume of complementary gabbroic cumulates roughly equal to the GRB erupted volume (150,100 km³, Kasbohm and Schoene, 2018) is implied. Detailed description of the mass balance calculation involved in this estimate is found in the Methods section.

This petrologically derived volume can then be compared to the estimated physical extent of intrusive CRBG volumes from mapped dikes encompassed by three named swarms: the Steens Mountain dike swarm (Steens Basalts), the Monument dike swarm (Picture Gorge Basalts), and the largest concentration of CRBG dikes: the Chief Joseph dike swarm (Imnaha; Grande Ronde, Wanapum, and Saddle Mountain Basalts; Fig. 1; Taubeneck, 1970; Fruchter and Baldwin, 1975; Tolan et al., 1989; Camp and Ross, 2004). This swarm was first recognized by Lindgren (1901) and described in detail by Taubeneck (1970). Dikes have been correlated to individual flows for some of the Wanapum and Saddle Mountain Basalts (Reidel and Tolan, 2013; Reidel et al., 2016). However, with a few exceptions (e.g., Petcovic and Dufek, 2005), physical correlations are more rare for the GRBs, in part because of the vast numbers of dikes and voluminous surface flows.

In this paper, we present an extensive CJDS data set based on notes and field maps compiled by Dr. William H. Taubeneck (WHT). This data set includes 4279 dike segments mapped at 1:24,000 scale across northeastern Oregon, southeastern Washington, and west-central Idaho; these segments define the structure of the CJDS in a regional context. The distributions of thickness, length, and orientation of dike segments help inform the following questions: What controls the orientation of dikes within the swarm? What are the implications for regional magmatic-tectonic relationships?

Our objectives are to: (1) describe the structural, textural, and chemical data set of WHT; (2) describe new observations to confirm and extend Taubeneck’s measurements through field work, remote sensing, and re-analysis of geochemical samples; (3) quantify dike swarm morphology and geometry; (4) characterize sub-swarm size and systematic trends; and (5) measure dike widths and characterize host-rock and dike interactions at the outcrop-scale. We then make use of these observations to (1) assess relative emplacement timing and crustal dilation; (2) assess controls of regional stress versus local mechanical variability on dike emplacement; (3) detail potential perturbations to the local and regional strain field from dike emplacement; and (4) investigate potential geochemical variations present within the dike swarm.

The CJDS stretches from the western Snake River Plain and northern Great Basin to the Columbia Basin. Most dikes are confined to an area ∼100 km wide by 350 km long (Figs. 1 and 2A; Supplemental Plate S1¹). The swarm traverses several lithospheric boundaries between Mesozoic accreted terranes and the margin of cratonic North America (Fig. 1). These terranes docked with North America by the Jurassic and were later intruded and sutured together by granitoids of the Wallowa and Elkhorn Mountains (Fig. 2A; Dorsey and LaMaskin, 2008; Schwartz et al., 2014a). However, dikes intrude through rocks that range from Proterozoic to Miocene in age.

During the Eocene, four notable events took place in this region of eastern Oregon. First, the Siletzia oceanic plateau accreted to the Pacific Northwest (Snavely et al., 1968; Irving, 1979; Wells et al., 2014). Second, Eastern Oregon experienced a period of exhumation, exposing granitoid bodies at the surface (Michaels et al., 2017). This unroofing in eastern Oregon was concomitant with exhumation of rocks within the Idaho Batholith to shallow crustal levels (Fig. 1; Giorgis et al., 2008; Fayon et al., 2017). Third, the 53–43 Ma Challis volcanics erupted just to the east (Gaschnig et al., 2016). Volcanism in eastern Oregon at this time period is linked with the dynamics of Farallon slab subduction (Liu and Stegman, 2011; Seligman et al., 2014; Camp et al., 2017b). The fourth event began with the eruption of the calc-alkaline Clarno units and early bimodal John Day Formation in central to eastern Oregon (McClaughry et al., 2009). Throughout this period of time, Oregon experienced south and westward increasing clockwise rotation and intermittent volcanism to the west in the ancestral Cascades arc (Wells et al., 2014; Camp et al., 2017b).

Volcanism in the Miocene began with the eruption of the Picture Gorge basalts at 17.23 ± 0.04 Ma (Cahoon et al., 2020), then by the Steens Basalts, starting ca. 16.9–16.6 Ma (Kasbohm and Schoene, 2018; Moore et al., 2018). There is also field and geochronologic evidence that after the eruption of the Steens Basalts, a southward progression in mafic volcanism into the northern Nevada rift began (Fig. 1; Camp and Ross, 2004).

The first erupted unit of the CRBG mapped in northeastern Oregon is the Imnaha Basalt, erupting at ca. 16.6 Ma contemporaneous with the last of the Steens eruptions to the south (Kasbohm and Schoene, 2018). This unit interfingers at its southernmost extent with the Steens Basalts (Hooper et al., 2002; Camp et al., 2017a). The Imnaha basalts have a volume of ∼11,000 km³ (∼6% of total CRBG volume) and erupted from vents in eastern Oregon (Taubeneck, 1970; Camp, 2013; Reidel et al., 2013; Kasbohm and Schoene, 2018).

The next oldest sequence of flows associated with the CRBG is the Grande Ronde Basalt (GRB). The GRBs are grouped into four paleomagnetic polarities: R1, N1, R2, and N2, the last three of which erupted in ∼70,000 years, with the entire Grande Ronde erupting between ca. 16.5 and ca. 16.1 Ma (Kasbohm and Schoene, 2018). These flows make up 72% of the total volume of the CRBG (Reidel and Tolan, 2013). Geochemically, they are classified as mostly basaltic andesites with 52%–57% silica versus 47%–53% silica present in Imnaha basalts (Hooper, 2000). Radiogenic isotopes suggest potentially significant assimilation of Idaho Batholith into the Grande Ronde basalts, a signal only weakly present in the Imnaha Basalts (Wolff et al., 2008). The bulk of the dikes within the CJDS are likely associated with the GRBs (Taubeneck, 1970).

Our database includes nearly 53 years of data, digitized from the collection of William H. Taubeneck, augmented by our own data collection during 2016, 2017, and 2018. To contextualize the CJDS, we also synthesized data on other dike swarms associated with the CRBG. Those data are detailed in the Discussion section.

Original maps from William H. Taubeneck (WHT) were scanned and then processed through ESRI ArcMap by georeferencing each 1:24,000 U.S. Geological Survey (USGS) quadrangle onto which WHT had mapped dikes. We digitized each linear, nonlinear, or curved feature marked on WHT’s maps, creating a shapefile containing all potential mapped dikes. Dike orientation, length, X and Y coordinates for each dike midpoint, and midpoint elevation were then extracted from the 4279 dike segments across northeastern Oregon, southeastern Washington, and west-central Idaho.

State geologic maps for Idaho, Oregon, and Washington were incorporated at this stage. The host rock was identified for each individual dike segment and recorded in the shapefile attribute table (Huntting et al., 1961; Bond et al., 1978; Walker and MacLeod, 1991). Host-rock lithologies were simplified into several categories: metasedimentary units (slates and marbles); metavolcanic units (greenstones and meta-rhyolites); granitoids (granite, tonalite, gabbro, and granodiorites); basalt (Miocene basalts or tholeiites); sedimentary (conglomerates, sandstones, and lake deposits), and andesite (andesite and intermediate volcanic rocks).

We then cross checked each dike segment location using 1 m resolution U.S. Department of Agriculture National Agriculture Imagery Program (NAIP) aerial photography. This allowed us to remove duplicate segments and those segments that are clearly not dikes, or non-CRBG dikes (older structures related to batholith intrusion in the Wallowa Mountains; Žak et al., 2012). These non-CRBG dikes are often light colored in the field and are visible in some satellite images and aerial photography. The final dike data set is available in the supplemental material for this publication as an ArcGIS shapefile².

Notebooks contained in the WHT collection document field observations associated with the mapped dikes in our database. Measurements in the WHT notebooks include: strike (n = 1891 measurements), dip (n = 167), dike thickness (n = 3289), elevation when the observation was made (n = 3043), date of observation (23 June 1953–2 September 2007), and approximate geographic location information (name of nearby peak or stream). WHT measurements of width likely represent a single locale along the length of a dike segment, so, at best, they should be treated as a “representative” width for the dike segment. WHT also included semiquantitative descriptions of dike and host-rock textures, such as vesicularity (n = 131), grain size (n = 194), inclusions (n = 117), weathering of dike (n = 37), partial melting of host rock (n = 120), and composition of dike (Grande Ronde versus Imnaha; Swanson et al., 1979). Each entry within the notebook that represented a dike was entered into an Excel spreadsheet. We quantified descriptions of dike vesicularity, inclusion density, and grain size on a 1–5 rating scale for each dike. This is an interpretation of Taubeneck’s notes. For example, dikes with a vesicularity of 1 contain no vesicles, while dikes with a rating of 5 correspond to “highly vesiculated” within Taubeneck’s notes. This rating scheme was also applied to the presence of dike-marginal melted wall rock, inclusions (e.g., xenoliths or other non-basalt lithologies), and dike grain size.

Inconsistencies in the data appear frequently. In particular, location data were difficult to follow. Taubeneck rarely lists township and range or other identifying coordinates. The best geographical information provided was elevation and general named locations (e.g., east of Glacier Lake; west and southwest of Eagle Cap; on west side of Hurricane Creek). Taubeneck’s system for locations on his maps (Loc. 1, 2, etc.) occasionally denoted multiple locations along a single dike, or one location for multiple dikes.

We have therefore been unable to clearly connect mapped dike segments directly with specific notebook entries. The digital database constructed from WHT’s work is included as a series of files in supplemental material for this paper (^{footnote 2}). The supplemental files also contain a table with dike midpoints, host rock, bearing, and length (Table S1³). Observations digitized from WHT’s notebooks and our own original observations made at CJDS dikes are contained in Table S2. Future work may provide a more robust link between these two data sets and extend the scope of dike mapping.

As part of an effort to check the accuracy of WHT maps to ensure any conclusions drawn from our database were robust, we have generated a new set of data to complement that of WHT. Our field observations and remote sensing aim to constrain relationships between fracture networks and dike orientation and the 3D geometry of dikes, both of which are poorly detailed in the WHT data set. Although our preliminary observations do not cover the entire CJDS area, they confirm the accuracy of the WHT data and contribute to a cohesive picture of multiscale structure of the CJDS described in more detail below.

Targeted field work took place in the Wallowa Mountains and a few limited areas to the south during the summers of 2016, 2017, and 2018. We found that WHT maps are, for the most part, easily verified in the field, and some also in remote sensing imagery. We suspect that WHT documented only the largest dikes, because many examples of smaller dikes seen in the field are not included in the WHT materials. We view extending the mapping and characterization of dikes as a promising direction for future work. The nature of dike tips and differences in dike–host-rock interaction are not directly addressed by this work.

We measured 768 joints in granitoid host rocks in the central Wallowas. Care was taken to measure joints greater than 10 m from any dike margin in an effort to avoid joints induced by dike emplacement (Hoek, 1991; Rubin, 1995). Dike width was measured with a measuring tape. Dike orientations were often difficult to measure on an outcrop scale due to the irregular nature of the dike–host-rock contact. To better constrain dike strike and dip in the field, the orientation of cooling joints within the dike that were normal to the dike margin was measured. The poles to these joint measurements were plotted on a stereonet and a best-fitting plane fit through those points, using the stereonet software (Cardozo and Allmendinger, 2013). This best-fitting plane provides a measure of dike orientation. This method was only applied to dikes that did not contain clear evidence of re-intrusion (e.g., selvage zones; Reidel et al., 2013) and had margin-normal columnar joints. The orientations of 40 dikes measured with this method are combined with noted strike-and-dip data from the WHT collection and orientations gathered through remote sensing.

Several remote sensing tools were used to provide data on dikes in areas that we did not visit, in addition to their use in verifying WHT mapped locations. We documented examples of: (1) crosscutting relationships; (2) dike strike and dip; (3) whether or not a dike was visibly segmented (e.g., en echelon behavior; Hoek, 1991). Crosscutting relationships between different dikes were observed in both Google Earth and 1 m resolution NAIP imagery.

Dike strike and dip were measured using 1 m horizontal resolution NAIP imagery in concert with Oregon Department of Gem and Mineral Resource light detection and ranging (LIDAR) (where available) and the 10 m resolution USGS National Elevation Dataset. A three- (or more) point planar-fit was performed using the LayerTools extension to ArcGIS (Kneissl et al., 2011). Dike segmentation was measured using Google Earth.

CJDS structures are scale-dependent. On the largest scales (>100 km), we compute statistics of the entire dike segment data set and fit distributions of dike segment average widths and lengths to commonly applied distributions to assess the degree to which the CJDS is comparable to other dike swarms (Krumbholz et al., 2014). To calculate the dike density, we used the line-density tool in ArcGIS, with a kernel window of 5 km (Silverman, 1986). The dike line density measured in km of dike per square km can be converted into dike density, assuming an average dike segment length (0.4 km). We calculated mean dike density by excluding areas with very low (0.1 km⁻²) to zero dike density, to isolate areas in which dike segments are exposed.

On smaller scales (1–100 km), we identify clusters of similarly oriented dike segments (within 20° bins) through a nearest neighbor analysis, comparing nearest neighbor distances of segment midpoints to mean and maximum segment length derived from the whole data set to demonstrate closeness of nearby segments (Fig. 3). While not a true clustering analysis, this procedure effectively highlights sub-swarm patterning of dike segments as a function of orientation, which can be compared to dike density (Fig. 4) At the smallest scale (< = segment length), we examine CJDS structures in outcrop.

W.H. Taubeneck sent numerous dike samples to the GeoAnalytical Lab at Washington State University for major- and trace-element analysis by wavelength dispersive X-ray fluorescence (XRF). We powdered, re-fused, and reanalyzed a number of original fused rock–Li tetraborate discs (“beads” hereafter) from the Washington State University bead archive on a ThermoARL Advant XP sequential XRF spectrometer using the methods of Johnson et al. (1999) (Table S3 [^{footnote 3}]).

Unfortunately, most of WHT’s samples cannot be map-located, so we cannot explore relations between dike chemistry and physical characteristics such as geometry or orientation. Instead, we compare the chemical data to those of the CRBG as a whole, and assign samples to CRBG formations on a preliminary basis.

The GRB lavas display remarkable homogeneity in some incompatible trace elements, in particular those that are least affected by the crustal contamination that accompanied crystallization-differentiation, such as Nb and Zr; relative 1 standard deviation variations are ∼10% (data of Wolff et al., 2008). Niobium is effectively incompatible in the crystallizing assemblage of Imnaha and GRB (plagioclase + pyroxene ± olivine) and is used as a fractionation monitor. It is also depleted in rocks of the Atlanta Peraluminous Suite of the Idaho batholith (Gaschnig et al., 2011) used by Wolff and Ramos (2013) as contaminants in their energy-constrained assimilation–fractional crystallization simulations (Bohrson and Spera, 2003) of the evolution from Imnaha to GRB compositions. Hence, fractionation can be tracked using Nb, and assimilation is ignored.

The Imnaha lavas are themselves somewhat evolved at 4%–7% MgO. Restoration of average parental Imnaha Basalt to 8% MgO, taken as the composition of basaltic liquid arriving in the crustal storage system from the mantle, using “reverse fractionation” gives a Nb content of 5.5 ppm. Reverse fractionation was done using MELTS (Ghiorso and Sack, 1995); average Imnaha Basalt composition was run at 10 °C below its liquidus at 1 GPa, f(O₂) = QFM (quartz-fayalite-magnetite), and the resultant assemblage (Ca-cpx) added back into the liquid in 10% increments until 8% MgO was achieved at ∼22% Ca-cpx addition. At 0.5 GPa, the MgO content of the pyroxene is little different, but plagioclase joins it on the liquidus; hence, the 1 GPa estimate is more conservative because less crystal addition is needed to achieve 8% MgO. Therefore, the calculated Nb in the starting magma is maximized and thus minimizes the volume of cumulates calculated below.

Data on dike swarm structures are presented below at telescoping scales from the order of the entire dike swarm (∼10s–100s km) to single dike segment scale (10s–100s m).

Dikes mapped at 1:24,000 scale by WHT within the CJDS are often dismembered by valleys, buried by younger deposits (fluvial or glacial), or obscured by active hillslope processes. As a result, our data set contains well-constrained dike segments but not total dike lengths.

Individual dike segments within the CJDS range in length from 30 m to 4.3 km. The distribution of dike segment lengths is well approximated as a log-normal distribution (Fig. 5A). The median length of a dike segment is 325 m. The log-normal mean and variance is 408 ± 8 m (root mean square error [RMSE] of 1.3*10⁻⁶), consistent across host-rock type (Fig. S2⁴). The median length for dikes hosted in granitoid rocks (n = 2291) is 363 m; for basalt 348 m (n = 1003); for metasedimentary rocks 243 m (n = 638); and for metavolcanic rocks 226 m (n = 346; Fig. S2). At the outcrop scale, there are some shorter dikes (<10 m in length) that branch off from a main dike. These probably are a small portion of the observed dikes within the swarm, but this is not quantified further here.

Dike width, as measured in the field, may not correspond to the active width of the dike during emplacement, which depends on dynamic pressure of magma flow. Measured width is usually assumed to correspond with one intrusion event, although dikes exhibiting evidence of multiple injections (from either multiple columnar joint sets or glassy internal contacts or selvage zones) are not uncommon within the CJDS (Taubeneck, 1970).

Widths vary somewhat along strike of many dike segments within the CJDS and appear to be somewhat depth dependent. Such variation is seen in other dike swarms and is expected near the dike tip (Jolly and Sanderson, 1995; Babiker and Gudmundsson, 2004; Paquet et al., 2007). In Figure 5B, we plot all Taubeneck dike widths (Table S4 [^{footnote 3}]). The median dike width across the CJDS is 8 m—comparable to other LIPs but much larger than typical dike widths in other settings (Babiker and Gudmundsson, 2004; Krumbholz et al., 2014). The minimum dike width from WHT notebooks is 0.2 m, although we have observed thinner dikes that are not included in the WHT data set. The maximum dike width from WHT is 146.3 m; however, no dike was observed during our field work approaching this size. The larger dikes are likely composite structures, recording multiple injections and selvage zones (Taubeneck, 1970; Reidel et al., 2013).

The distribution of dike widths has been proposed to inform processes associated with emplacement; for example: (1) magmatic overpressure—Babiker and Gudmundsson, 2004; (2) length of time for flow within the dike—Parfitt and Wilson, 1994; (3) depth within the magmatic system—Delaney and Gartner, 1997; Woods et al., 2006; and (4) host-rock rheology and elastic moduli—Krumbholz et al., 2014; Karlstrom et al., 2017. A Weibull and a log-normal distribution both fit the swarm-scale distribution of dike widths well (Fig. 5B inset). This has been found in other settings (Jolly and Sanderson, 1995; Delaney and Gartner, 1997; Krumbholz et al., 2014).

The mean dike line density excluding the lowest values (<0.1) is 0.5 ± 0.37 km of dike per square km. This corresponds to a mean dike density of ∼1.2 ± 0.9 dikes per square km. Maxima in dike density occur in the Wallowa Mountains and several of the intrusions to the south, reaching 1.5–1.9 km of dike per square kilometer (Fig. 4). Assuming an average dike length of 0.4 km, this maximum dike density results in five dike segments per square kilometer. Dike density is also greatest in areas that also have the greatest variability in dike orientation (Fig. 4A).

There are several identifiable trends in dike orientation across the CJDS. The dominant trend is NNW (median bearing of −352°; Fig. 2D). This trend is generally consistent across all lithologies in which dikes are observed (Fig. S1 [^{footnote 4}]). The median bearing for dikes in granitoids is 351°; in metasedimentary rocks is 355°; in metavolcanic rocks is 351°; and in basalt is 351° (Fig. S1). The dominance of this orientation makes it difficult to perceive other trends that appear visible in map-view (Fig. 2A; Plate S1 [^{footnote 1}]).

A moving average rose diagram with a 6° window size (Fig. 2D) suggests smaller, yet still significant, trends within the dike swarm. This is confirmed by nearest neighbor analysis (Fig. 3A), which reveals that segments at all orientations are commonly observed with similarly oriented segments nearby. Clusters of dike segments within the primary dike orientation (NNW) dominate: on average, segments within 10 degrees of N-S exhibit ∼70 nearby segment midpoints (maximum of 211) closer to the CJDS maximum segment length of 4.3 km (Fig. 3A). In contrast, dikes within 10 degrees of EW have on average 11 segments with midpoints this far away (maximum 27). Remarkably, this trend still holds when the nearest neighbor threshold is reduced to 0.4 km (the mean CJDS segment length), demonstrating that some areas of the CJDS exhibit extremely dense clusters of segments.

Figure 3B shows one such region in the central Wallowa Mountains; this region exhibits dense and overlapping dike segment clusters of different orientations. Dikes are colored by their orientation broken into 20° angular bins. The observation that dike segments are often proximal to other segments with a similar orientation suggests that they may in some cases reflect a common deeper structure, broken up in map view by topography or covered by younger deposits.

Dike dips were collected from the WHT database (n = 112). These data were combined with original field measurements (n = 38) and remote sensing measurements (n = 172). There is a large range of dips. The lowest angle intrusion (sill) reported dips at 4° (from WHT). The steepest dikes are vertical (Fig. 6). There is significant variance within dike dips but most are subvertical; the overall distribution of dips is negatively skewed with a median dip of 65° (Fig. 6). Dikes also appear to dip at different orientations depending on the host rock. Those hosted by granitoid rocks appear steeper, with a median dip of 68°. Dikes hosted in metasedimentary rocks are less steeply dipping, with a median of 55°. Of the total dike dip measurements, data were available to classify 174 dikes as east-dipping and 128 dikes as west-dipping.

Only 17% of all intrusions measured within the CJDS (n = 322) have dip less than 45° (Fig. 6). This is different from some other LIPs. For instance, the bulk of exposed intrusions within the Siberian Traps are sills, which make up 16% of the total areal extent of the Traps (∼10⁷ km²; Svensen et al., 2009; Burgess et al., 2017). Sills up to 22 m thick are also found within the Deccan Traps, another LIP where dikes appear as the apparent dominant intrusive style (Sheth et al., 2009; Duraiswami and Shaikh, 2013).

The WHT samples span beyond the distribution that generally defines the range of CRBG lavas (Fig. 7), but most samples fall within that range, so in principle can be assigned to a unit using the chemical correlation approach of Hooper (2000). Of 93 samples total, most have affinities with Grande Ronde Basalt (Fig. 7). There is some overlap between Grande Ronde and Imnaha Basalt in the 52%–53% SiO₂ range (Wolff et al., 2008; Wolff and Ramos, 2013) with current methods. Five andesites are more silicic than any known CRBG lava (59.1%–61.2% SiO₂) but have incompatible elements consistently elevated by ∼25% over typical Grande Ronde lavas (Table S3 [^{footnote 3}]), hence could plausibly be related differentiates.

Of uncertain affinity is a group of four samples (WT28-79, WT28-847, and WT28-913) that are more magnesian (∼9% MgO) than any previously analyzed CRBG samples with the exception of a few Steens Basalt flows (Johnson et al., 1998; Moore et al., 2018) and Monument dikes (Bailey, 1989). No lava that erupted from the CJDS with such high MgO contents has been found, although the basalt of Robinette Mountain (Eckler Mountain Member, Wanapum Basalt) is generally similar (see Table S3 and Hooper, 2000).

The average Nb content of the GRB is 11.4 ppm (data of Wolff et al., 2008), corresponding to 52% crystallization of a parental Imnaha Basalt. Using volume-weighted average compositions of each of the four magnetic polarity divisions of the GRB (Reidel and Tolan, 2013) gives the same result. Without taking the volumetric contribution of crustal contamination into account, this yields a volume of gabbro added to the crust of ∼160,000 km³. Wolff and Ramos (2013) estimate the recycled crustal contribution to the GRB at 9%–59%, yielding net cumulate additions to the crust of 147,000 km³ and 66,000 km³ gabbro, respectively, with a mean of 107,000 km³. This volume has not been compared to estimates of the intrusive volume of the CRBG. Herein, we endeavor to do so.

Dike segmentation has a variety of origins. It may occur during emplacement due to rotations of remote stresses along the dike propagation path, called en-echelon segmentation (Fig. 8B; Pollard et al., 1982). Dikes may also exhibit branching or break-out structures that segment observed outcrops (Fig. 8A; Hoek, 1991). Inelastic deformation either syn- or post-emplacement may also segment dikes (Mathieu et al., 2008).

Dikes within the CJDS exhibit multiple apparent segmentation styles. There are many segments of dikes that are approximately linear over their entire length. However, a small percentage (<5%) of all dikes that we observed via satellite imagery contain segments that appear en echelon or “zig-zag” as classified by Hoek (1991) (see our Fig. 8B). There is no clear preference between left-stepping (n = 91) or right-stepping (n = 104) segments, over the dikes that we have analyzed (Fig. S3 [^{footnote 4}]). Locations of mapped en echelon offsets are viewable in the Supplemental Google Earth File⁵.

A number of the CJDS dikes were continuously active for an extended period of time as they fed surface CRBG eruptions (perhaps several years; e.g., Swanson et al., 1975; Thordarson and Self, 1998; Petcovic and Grunder, 2003; Petcovic and Dufek, 2005; Karlstrom et al., 2019). Many other dikes were discontinuously active and reoccupied by multiple magma injection events, as seen by internal selvage zones or multiple internal chilled margins (Taubeneck, 1970; Reidel et al., 2013). We use crosscutting relationships to document the progression of diking within the CJDS. We measured 46 crosscutting relationships between dikes exposed in the central Wallowa Mountains and Cornucopia region through a combination of field and remote-sensing observations to generate an aerially extensive data set (Fig. 9; Table 1; Google Earth File [^{footnote 5}]). Observed crosscutting relationships suggest that the NNW-oriented dike segments are the oldest. However, in the region analyzed, these do not constrain the relative ages of subsequent diking, which may have been oriented NE-SW, NW-SE, E-W, or N-S.

The WHT notebooks contain 3113 entries with dike elevation and dike thickness noted (see Table S2 [^{footnote 3}]). Of this data set, 1731 measurements occur above 1500 m. The bulk of the topography in northeastern Oregon above 1500 m is exposed in the Wallowa Mountains, which are known to have uplifted post-CRBG and provide a depth-section through the magmatic plumbing system (Hales et al., 2005; Schoettle-Greene and Duvall, 2016). There are Imnaha basalt flows atop ∼3000 m peaks with ∼1500 m of relief below these summits (Hales et al., 2005). So, while we cannot robustly identify the geographical location of WHT width measurements, we can test a relationship between dike elevation (as a proxy for paleodepth) and thickness in the Wallowas.

A linear regression through available data on dike elevation suggests that dikes thin toward shallower paleodepths in the magmatic system (regression slope of −0.002; Fig. 10). The adjusted R-square reveals a poor fit between elevation and dike thickness (Adj. R² = 0.02 explains only 2% of the variance in dike width; Zar, 2010). However, a P-value for this fit indicates the slope has significance (P-value = 4*10⁻¹⁰; Fisher, 1992; Zar, 2010).

Taubeneck (1990) hypothesized that dikes within the CJDS interact with and perhaps follow preexisting host-rock joints. Dikes in other areas have been observed both to correlate with host-rock joints (e.g., Jolly and Sanderson, 1995; Delaney and Gartner, 1997) or to propagate without regard for host-rock joints (e.g., Delaney et al., 1986). To further evaluate the importance of preexisting host-rock fractures in the CJDS, the orientations of 251 joints were measured in the central Wallowa Mountains. Care was taken to measure joints >10 m from any dike to ensure they did not form as a result of diking. Some dikes hosted within granitoid rocks appear to have some relationship with the preexisting host-rock fractures on a <10–100 m scale (Fig. 8A). Either the dikes are following and exploiting these preexisting fractures (e.g., Delaney and Gartner, 1997), or the fractures are forming as a result of diking (Delaney et al., 1986).

These joint orientations are compared to the orientation of 185 dike segments exposed in the same area (Fig. 11). There is an observed correlation between the 340° dike orientation and joints with a similar orientation, and the 60°–80°–oriented joints correspond with the 60°-oriented dikes (Fig. 12B). This observation is consistent with Taubeneck (1990), suggesting that some dikes interact with host-rock joints—at least locally.

Because these joints were observed >10 m from dike margins, dikes were not forming the joints but likely intruding preexisting fractures. It is rare to directly observe dike material intruding host-rock joints in the CJDS (both in WHT notes and our observations). The dike exposed in the north face of Eagle Cap Peak appears to be an exception with a break-out segment that follows the orientation of host-rock fractures (Fig. 8A). However, dikes that do not follow fractures also occasionally occur (Fig. 11B), such as flame-like fracture sets suggestive of hydrofracture (see an example from the Cornucopia region of the Wallowas in Fig. S5 [^{footnote 4}]). At a larger 1–10s of km scale, CJDS dike segments and segment clusters cut across fractures, host-rock changes, and even crustal boundaries between accreted terranes (Fig. 2A) and still appear linear.

The WHT database and adjoining data compiled here provide a detailed snapshot of the shallow CRBG plumbing system, a cross section through a massive continental dike swarm that is unique in preservation, and exposure. Here, we synthesize these data to generate a magma-tectonic framework for understanding the CJDS. In keeping with the structure of the Results section, we move progressively inward from province to outcrop scale. We first compare the CJDS to other mapped dikes within the CRBG, then discuss swarm-scale trends in segment geometry. We estimate the total intrusive volume of the CJDS and discuss implications for crustal magma reservoirs associated with the CRBG. Dike swarm structures are used to infer changes through time of the principal stresses in the shallow crust, regional dilation associated with dike emplacement, and possible swarm-scale dike thinning toward the surface. Finally, we connect segment geometry with simple models for magma flow to argue that most observed CJDS dikes were not feeders to surface flows.

While the CJDS is the largest dike swarm in both aerial extent and number of dike segments, three other CRBG swarms have also been identified throughout eastern Oregon, southeast Washington, and west-central Idaho: the dikes that erupted the Steens Basalts; the Monument dike sub-swarm that erupted the Picture Gorge Basalts, and the Ice Harbor dike sub-swarm that fed the Ice Harbor member of the Saddle Mountain Basalts (Fruchter and Baldwin, 1975; Reidel et al., 2013; Blakely et al., 2014; Camp et al., 2017a; Kasbohm and Schoene, 2018). We have compiled a map of all known CRBG dikes here to contextualize the CJDS, and we detail below where these data are sourced from (Plate S1 [^{footnote 1}]; Table S5–S7 [^{footnote 3}]). We also created a simplified map of the dikes related to the CRBG, for use in plotting an accurate representation of dike locations, by removing excess vertices along dike segments and randomly down-sampling the population of dikes >1 km in length (Plates S2 and S3).

The Steens Basalts dikes are exposed on the eastern flank of Steens Mountain. We mapped 69 basaltic dikes using Google Earth Imagery. These dikes are highlighted in Plate S1 (^{footnote 1}) and are available as a shapefile (^{footnote 2}) in the supplemental material. We also incorporated dikes from the Picture Gorge Basalts (Monument dike sub-swarm) exposed in central Oregon. These dikes were mapped originally by Brown and Thayer (1966) and Wallace and Calkins (1956), and then digitized by Cahoon et al. (2020). The Monument dikes are also available as a shapefile and are plotted in Plate S1.

The youngest set of dikes mapped as part of the CRBG magmatic event are those inferred to be from the Ice Harbor flows (8.8–8.5 Ma; Mann and Meyer, 1993; Hutter, 1997; Reidel et al., 2013). There are some limited exposures of Ice Harbor dike sub-swarm in the Pasco Area (Fig. 1; Plate S1; Hutter, 1997), but for the most part, these ∼130 dike segments were imaged using high-resolution aeromagnetic data (Plate S1; Blakely et al., 2014). We digitized the inferred dike positions and integrated this data set into our collection of CRB-related intrusions. These data are available in a shapefile (^{footnote 2}) and visible in Plate S1 (^{footnote 1}). The Ice Harbor Dikes are the only features within this data set that have an inferred position; all other dikes are constrained by outcrops.

Length distributions of these other dike sub-swarms provide context for the CJDS. Steens dike segments have a median length of 189 m, with a maximum length of 1.7 km and a minimum length of 26 m. Monument dike segments have a median length of 1055 m, with a maximum length of 3.8 km and a minimum length of 260 m. Ice Harbor dike segments have a median length of ∼3400 m, with a maximum length of 22 km and a minimum length of 593 m (Fig. 5B; Plate S4). That the Ice Harbor dikes are observed to be significantly longer than the other swarms may be explained by their observation through geophysical means. The Ice Harbor segments are still largely subterranean and have not been dissected by topography or obscured beneath younger deposits as is the case with the other dike segments in the CRBG. In fact, it might be appropriate to consider these dikes as a model for the “primary” CRBG dike segment length distribution before erosional dissection and sedimentary cover. We note that these dikes are from the Ice Harbor Member, which is relatively small in volume compared to other CRBG formations (75 km³; Reidel et al., 2013).

Our database across thousands of dike segments, transgressing different episodes of the CRBG event, also allows for a regional assessment of dike orientation (Plate S1 [^{footnote 1}]). The dike swarms that represent the oldest portions of this magmatic event are the Steens (16.6–16.5 Ma), Monument (17.2–16.06 Ma), and perhaps some number of the NNW-oriented dikes of the CJDS, which are crosscut by many other dike orientations (Fig. 9; Kasbohm and Schoene, 2018; Cahoon et al., 2020). Dikes of the Monument swarm are dominantly oriented NW-SE (330° ± 14°); Steens dikes are dominantly oriented NNE-SSW (10° ± 27°), and the oldest trend identified from our preliminary field work in the CJDS are the NNW-oriented dikes (Fig. 2; Plate S1). The youngest dikes mapped in our database, those associated with the ca. 8.5 Ma Ice Harbor flows, are dominantly oriented NW (336° ± 12°). This is a notable ∼10°–15° shift westward from the dominant trend of the CJDS, which is 352°.

Previous authors identified a radial orientation across these different dike swarms, invoking a rising plume head to explain dike orientation and dike distribution (Ernst and Buchan, 1997; Glen and Ponce, 2002; Camp and Ross, 2004). Such a model would require large-scale lateral transport of magma through the dike swarm (e.g., Ernst and Baragar, 1992), or alternatively, migration of magmatism could be accomplished through lateral migration of localized magma chambers (e.g., Camp and Ross, 2004; Sleep, 2008). However, a lack of palinspastic reconstructions complicates the hypothesis of radial diking in the case of the CRBG dikes. Extension and rotation have taken place since the Miocene in northern Nevada and southern Oregon (McQuarrie and Wernicke, 2005; Camp et al., 2015). Tectonic deformation is less well constrained in the areas of northeastern Oregon and southeastern Washington, where most of the CJDS is exposed. Work conducted by McCaffrey et al. (2013) and Wells et al. (2014) places poles of rotation for post-Miocene extension and rotation in northeastern Oregon. This pole of rotation implies little tectonic extension has taken place in northeastern Oregon since the Miocene; the lack of strong tectonic gradients in the region of the CJDS supports a magmatic stress–dominated system.

Nearly 50% of all dikes exposed across the CJDS are hosted within granitoid bodies. Taubeneck (1970, 1997) suggests that dikes preferentially intruded this rock type. To test WHT’s hypothesis, we examined the hypsometry of different rock types across northeastern Oregon, comparing these elevations with the elevations of exposed dikes (Fig. 12). The median dike elevation is ∼1900 m; the median elevation in northeastern Oregon is ∼1200 m (Figs. 3B, 12A, and 12B). The two most common host rocks—basalt and granitoids (Figs. 12C and 12D)—have median elevations of ∼1200 and ∼1750 m, respectively. The distribution of dikes below ∼1900 m is constant with elevation; however, there is a notable increase in the occurrence of dikes above 1900 m. This elevation is closer to the elevation of the Pleistocene equilibrium line altitude (ELA; ∼2000 m) than any other topographic and lithologic metric (Fig. 12A; Meyer et al., 2004).

If dikes were focused toward granites, then one would expect to see a peak in dike occurrence at similar elevation areas dominated by granitoid rocks (1500–2250 m). However, dikes still occur with increasing probability at higher elevations (Fig. 12A). We posit that more dikes are exposed at higher elevations due to glacial erosion, rather than due to a true increase in dike abundance in granitoids. A full comparison between dike elevation, hypsometry, and all mapped host-rock types is shown in Figure S6 (^{footnote 4}). The exposure bias of dikes in higher elevation areas also suggests that there may be numerous unexposed dikes at lower elevations and similar paleodepths.

Taubeneck (1970) extrapolated a maximum number of dikes from the highest dike density within the Wallowas across the extent of the CJDS. He suggested that as many as 21,000 dike segments may exist across the region, with a potential upper limit of 30,000 segments. With a more complete database than was available in 1970, we can replicate this calculation. Within the main area of the CJDS (100 × 350 km), the greatest dike density measured within the Wallowas (approximately five dike segments per km²) yields ∼175,000 dike segments across this region, well beyond WHT’s estimate (Fig. 4B). Because the Wallowas represent ∼10% of the exposed area of the CJDS, this is likely a significant overestimate. However, extrapolating from the mean value of dike density (∼1.2 ± 0.9 dikes per km²) gives ∼42,000 dikes across the CJDS, which is also larger than WHT’s estimate. Of course, because dikes are clustered, assuming homogeneous dike density across the swarm is an upper bound on the number of dike segments.

Dike intrusion represents unrecovered elastic strain normal to the propagation direction, with smaller contributions from thermal and mechanical erosion of wall rocks (Anderson, 1942; Hoek, 1991; Parfitt and Wilson, 1994; Fialko and Rubin, 1999). Integrated across a large dike swarm such as the CJDS, such opening implies significant dilation of the upper crust during emplacement of the CRBG.

In Equation 1, D_k is the component of dilation (in meters) oriented in the k = N-S, E-W direction for the i th bin of either longitude (for N-S dilation) or latitude (for E-W dilation). y_i is the start of the longitude and latitude bin, y_i+1 is the end. We chose a bin width of 1.7 km in each case and tested that the results are not highly sensitive to bin size. is one of M total dike segments within each bin. It is scaled to a common origin, dotted into either N-S or E-W unit vectors ⁠, then multiplied by median dike thickness μ to obtain N-S or E-W dilation in meters.

This procedure results in a profile of average dilation that is plotted in Figures 2B and 2C. This is a lower bound for crustal dilation in the CJDS, for several reasons: (1) surface exposures of many dikes are likely covered by CRBG lavas or post-Miocene sediments; (2) dike segment spatial density likely increases with depth; and (3) many smaller dikes are not well documented in the WHT data.

In the observed pattern of dilation (Figs. 2B and 2C), E-W extension peaks in the Wallowa Mountains (∼45.16°N, ∼117.30°W), the area of greatest dike exposure with ∼500 m of dilation. N-S dilation peaks in the same location with ∼150 m of dilation, reflecting the large dispersion in dike orientations here (Fig. 4A). The envelope of dilation decreases substantially away from this central location.

Just like estimating the total number of dikes, we can estimate the near-surface strain associated with the CJDS by scaling dilation over the 100 km E-W and 350 km N-S area over which much of the CJDS is exposed. Maximum strains associated with peaks in E-W and N-S dilation are ∼0.5% and ∼0.06%. A maximum amount of strain can be calculated if the greatest dike density observed in the Wallowa Mountains (approximately five dike segments per km²) is scaled to the 100 × 350 km area over which the CJDS is exposed. This calculation yields ∼13% E-W–oriented strain or 2.7 km of dilation, and ∼0.3% N-S–oriented strain or ∼1.3 km of dilation. To put this scale of dilation in perspective, John et al. (2000) calculated ∼1 km of extension due to diking across the contemporaneous northern Nevada rift. Increased precision of dating suggests that three magnetostratigraphic divisions of the Grande Ronde basalts erupted in an increasingly short window of time ∼70,000–400,000 years (Kasbohm and Schoene, 2018). Given that a large fraction of the CJDS dikes is Grande Ronde in composition, such rapid emplacement indicates that the transient EW dilation rates represented by the strain calculated here are comparable to modern plate speeds (∼cm yr⁻¹)!

We examined dilation in cardinal directions (E-W and N-S) in Figures 2B and 2C. However, the dominant dilation direction is slightly rotated from E-W (Fig. 2D). While the bulk of dike dilation appears consistent with other structures present in the intermountain west (i.e., the northern Nevada rift; Fig. 1; Colgan, 2013; Camp et al., 2015), CRBG-related crustal dilation oriented in other directions has not previously been studied (Figs. 2D, 3, and 13; Colgan, 2013). Although the relative contribution of strain in these directions is small, it indicates that either (1) dikes are not responding purely to an ENE-WSW–oriented minimum compressive stress direction, or (2) the stress field changed through time. It is also important to note here that our observations constrain the stresses within the upper ∼2 km of the crust. The stress field within the mid-crust and the orientation of deeper intrusions may be different (Isola et al., 2014; Corti et al., 2019).

In Equation 2, σ_t is the tensile strength of the rocks, generally in the range of ∼1–10 MPa. With Young’s Modulus E ∼30 GPa, if H ∼50–100 km to match the CJDS dimensions and u ∼ 1 cm yr⁻¹ to produce dilation of ∼1 km over ∼100 k.y., then I ∼10⁻²–10⁻⁴ yr⁻¹ from Equation 2. Dike intrusion at a higher rate, plausible given the larger number of segments in the high-density regions of the CJDS, would transiently wipe out an extensional stress field as is suggested by clusters of dike segments with overlapping orientations (Fig. 3A). Equation 2 balances extensional tectonic stresses with intrusion. If we estimate the extensional stress associated with dike intrusion alone, simply scaling the stress associated with dike intrusion: ~3 MPa, where 10⁻⁴ is the strain associated with a single dike. Given that the deviatoric stress in the upper crust ranges from 20 to 30 MPa (Stüwe, 2007), the intrusion of even 10–20 dikes in the absence of extension could alter regional principal stresses. The observation of a dominant early orientation for the CJDS segments, which is then crosscut by dike clusters with variable orientations (Figs. 2D and 9), suggests that dikes intruded in a predominant NNW direction until reaching an isotropic (e.g., Muirhead et al., 2014), or ortho-rhombic condition with unequal principal axes of stress (Reches, 1978).

If such a transition from tectonic to magmatically dominated stresses occurred, crustal magma transport could be increasingly affected by the CRBG-imposed magmatic reservoirs (e.g., Karlstrom et al., 2009, 2015) as has been hypothesized by Wolff et al. (2008) for the CRBG, or stresses associated with surface loads from the lava pile (Burgess et al., 2017). Dikes in this scenario could be more easily affected locally by time and spatially varying magmatic stresses, but re-orient farther afield to a tectonic stress regime (Muller and Pollard, 1977). This would also support the observed dikes following joints in one stress condition and perhaps at another time cutting their own path regardless of joint orientation (Fig. 11; Delaney et al., 1986). Such behavior might be used in the CJDS to map the spatial distribution of deep magma reservoirs, if vertical versus horizontal flow can be established (Fig. 4A; Muller and Pollard, 1977; Muller, 1986; Jolly and Sanderson, 1997). Moreover, the behavior of CRBG dikes cutting their own path is consistent with observations made by Colgan (2013) and Camp et al. (2015) in the northern Nevada rift, where diking is offset by ∼45° from the predominant Miocene extensional direction.

The high rates of deformation implicated by dike intrusion may be exacerbated by unique conditioning of the lower crust above a mantle plume head. Much work has been done to discuss the potential movement of lower-crust and lithospheric mantle near a spreading mantle plume head in observational studies (e.g., Camp and Hanan, 2008; Darold and Humphreys, 2013) and modeling studies (e.g., Burov et al., 2007; Guillou-Frottier et al., 2007; Cloetingh et al., 2013). If plume-head spreading resulted in lateral or vertical movement of the lithospheric mantle and lower crust beneath the CJDS, then dike-induced dilation of the upper crust may not be directly communicated to the lower crust (Lustrino, 2005).

Here, v is the Poisson’s ratio of the host rock, L is the length of the dike (perpendicular to width, and P_o is the magmatic overpressure (pressure in excess of lithostatic pressure) within the dike. Dike overpressure depends on a range of factors including magma buoyancy, magma reservoir pressure, and the differential background stress (Gudmundsson, 1990; Geshi et al., 2010).

Assuming that the median CJDS widths are representative of w in Equation 3, we can estimate the horizontal lengths expected in the CJDS. Using a range of magmatic overpressures (P_o = 1–10 MPa), with E = 30 GPa, and v = 0.25 mean dike lengths range from L = 10–100 km. This is similar in magnitude to the largest remotely sensed Ice Harbor dike length but larger than the median CJDS segment lengths by a factor of ∼100. An LEFM scaling thus supports our hypothesis that Ice Harbor dike segments are more representative of dike length than CJDS segments. Of course LEFM in its simplest form is likely not a complete description of LIP-scale dikes, and models that better account for multiple interacting hydraulic fractures (e.g., Bunger et al., 2013) should be applied.

The WHT data suggest a small but robust anticorrelation between dike width and paleodepth within the CJDS (Fig. 10). There are three possible explanations for this trend: (1) Near the dike tip region, the fracture will narrow (Geshi et al., 2010). Therefore smaller widths at shallower depths could be an observation of long-wavelength tipward dike narrowing. (2) Viscous pressure loss from a magma flowing away from a chamber within an elastic-walled dike results in thinning, although this does depend on the rate of viscous pressure drop and spatial distribution of stresses in the medium (Woods et al., 2006; Pinel et al., 2017). (3) Dike thinning is based on changes in buoyancy of the magma relative to the host rock. Where the density of the crust matches that of the magma in the propagating dike, the level of neutral buoyancy, the dike may either change direction and propagate laterally or stop propagating (perhaps temporarily) and form a magma reservoir (e.g., Sparks et al., 1984). The dike will be thickest at its level of neutral buoyancy and thin toward the surface (Lister and Kerr, 1991; Townsend et al., 2017).

One possible explanation for this apparent lack of sills within the CJDS is that a shallower paleodepth within the magmatic plumbing system is exposed than in other LIPs. In this case, future seismic work may detect sills or sill complexes at depth. Moreover, discontinuous bedded units in an extensional environment with a horizontal minimum principal stress do not favor sill formation, which could also explain the difference between the CRBG and other provinces (Menand et al., 2010). The most regularly bedded rock types in the region are the flows of the CRBG; older host rocks in the CJDS area have undergone significant postdepositional deformation (White et al., 1992; Dorsey and LaMaskin, 2008; Stanley et al., 2008; Žak et al., 2012, and references therein).

The extrusive component of the CRBG has been estimated as 210,000 km³ (Reidel et al., 2013). With the WHT data, we can now estimate the intrusive volume of the CRBG contained within the CJDS. Estimates for the ratio between intrusive and extrusive (I:E ratio) rocks in magmatic provinces range from 1:1 up to 16:1 (Shaw, 1980; Wadge, 1980; Crisp, 1984; White et al., 2006). More recent work in South America suggests the ratio may be as high as 35:1 (observed geophysically) or 75:1 (from petrochronology; Tierney et al., 2016; Ward et al., 2017). These estimates vary depending on tectonic setting and mantle melting regime, and could change over the duration of a magmatic system (Karlstrom et al., 2017; Colón et al., 2018). Below, we develop several potential models for the intrusive to extrusive ratio for the CRBG.

The source depth of CJDS dikes is unknown, which complicates an accurate estimate of dike volumes. Geophysical evidence points to cumulate layers at the base of the crust in the CRBG (Catchings and Mooney, 1988; Davenport et al., 2017), resulting in greater seismic velocities that are a common occurrence in flood basalt provinces worldwide (Ridley and Richards, 2010). The presence of large cumulate reservoirs is likely required for mantle-derived ultramafic melts to fractionate into basalts (Karlstrom and Richards, 2011). However, thermobarometry suggests CRBG magma crystallization throughout the crust rather than just at the Moho (Rodriguez and Sen, 2013). Moreover, significant crustal assimilation is inferred for the GRB formation (Wolff et al., 2008). Depth distribution of magma reservoirs for the CRBG thus remains somewhat of an open question.

We therefore treat two magma reservoir depth scenarios—a Moho-level dike source at 30 km and a mid-crustal dike source at 15 km—to calculate the total volume of the CJDS from observed dike segment densities. If we simply sum the total length of mapped dikes in our data set (∼2000 km length) and assume a width of ∼10 m, we find a total intruded volume of ∼250–450 km³. However incomplete exposure suggests this is a lower bound. If we instead take the mean mapped dike density (∼1 dike/km²) as representative of the majority of the CJDS (defined by our 100 × 350 km rectangle), this results in ∼35,000 dike segments. With a representative length per segment of ∼1 km, we get ∼5–10.5 × 10³ km³ intruded volume for the two depth scenarios. The highest estimate permitted by our data arises if we assume the maximum observed dike spatial density (approximately five dikes/km²) is a representative crustal average. This results in 3–6 × 10⁴ km³ for the CJDS intruded volume.

The heterogenous nature of dike density throughout the swarm likely makes these estimates upper bounds to CJDS volume. Indeed, two orders of magnitude difference in volume estimates for the CJDS reflect our considerable uncertainty both in spatial density across the swarm that is not exposed, and to a lesser extent, the unknown source depth. The range ∼2.5 × 10²–6 × 10⁴ km³ represents 0.1%–34% of the entire CRBG volume, the upper limits of which are surprisingly large given that we have neglected the presumably larger reservoir volumes that must have sourced the dikes. This volume of spatially distributed intrusive flux likely has a significant impact on crustal rheology (Karlstrom et al., 2017; Perry-Houts and Karlstrom, 2019) and contributes non-negligible volumes to the overall volatile budget and associated atmospheric impacts of the CRBG (Self, 2006; Black and Manga, 2017). For example, if we assume 0.1–0.2 wt% primary SO₂ concentrations as inferred from melt inclusions in Grande Ronde near vent tephra (Davis et al., 2017), this volume accounts for 1s–100s Gt SO₂ and is comparable to the largest CRGB flow units.

Our petrologic mass-balance calculations suggest a maximum of ∼160,000 km³ of gabbro were added to the crust during the eruption of the Grande Ronde Basalts (similar to the extruded volume of the GRBs). This volume of basalt, if added to the entire region covered by the CJDS (100 × 350 km), would contribute 4.5 km to the thickness of the crust. If we assume that storage and crustal assimilation occurred in a more localized region, the resulting crustal thickening increases correspondingly (Wolff et al., 2008; Wolff and Ramos, 2013). Current seismically derived crustal thickness maps for this region do not support such significant thickening (Davenport et al., 2017). Combining geophysical mapping of crustal structure with the dike structural evidence for crustal dilation associated with dikes seems a promising direction for reconciling this apparent discrepancy. This may in turn help refine the crustal extent of the CJDS.

Additionally, if the range of volumes calculated exclusively for the CJDS is compared to the volume for the corresponding surface flows, the intrusive to extrusive ratios are: 1:1000 to ∼7:20. These ratios are significantly lower than the ratios seen in other volcanic provinces (White et al., 2006; Tierney et al., 2016). When the petrologically derived intrusive volume of ∼160,000 km³ associated with the GRB is used, the intrusive to extrusive ratio is ∼1:1.

In shallower portions of the CJDS, some dikes have been robustly linked to surface flows through geochemical correlation (e.g., Reidel and Tolan, 2013; Reidel et al., 2016). In the central Wallowa Mountains, the deepest exposure of the CJDS, there are poor links between dikes and flows and very little published information on how long feeder dikes may have been active. Karlstrom et al. (2019), extending the petrographic study of Petcovic et al. (2001), used low-temperature thermochronology and thermal modeling to show two dikes within the CJDS have different thermal histories that imply variable flow durations. One ∼8-m-wide dike within the Wallowa Mountains in the CJDS was long lived (active flow duration perhaps as long as ∼10 years), whereas a similar dike to the south in the Cornucopia region actively transported magma for at most ∼2 years and likely less (Karlstrom et al., 2019). The WHT database allows us to coarsely extend these case studies, because only ∼3% of dikes in the database have some degree of partial melt at the dike margin that indicates long-lived transport. In fact, an inference based on simple flow models described below requires that most dikes were not feeders to the CRBG eruptions.

To match estimates of CRBG from geochronology, ≤400,000 years for eruption of the GRB, most dikes of the CJDS cannot be feeders (Kasbohm and Schoene, 2018). Moreover, it is likely that flow was localized to smaller regions of those active feeders, on the order of tens of meters, not kilometers (e.g., Wylie et al., 1999). Using the same parameters enumerated above yields active length L = 0.5 m to match a flux of F = 5000 m³ s⁻¹ roughly equivalent to the estimated eruptive flux from the 1783–1784 Laki eruption in Iceland (Thordarson and Self, 2003). An active width less than the dike width is also physically unlikely, unless a fraction of the observed width actively transported magma. This change in width would have significant effects on the calculation above of the overall time period for CRBG eruption. We view reconciling the observed geometry and extent of the CJDS with the observed erupted volume of the CRBG a major challenge for future work.

The WHT data set contains many qualitative descriptions of dike textures, inclusion density, grain size, and vesicularity, tabulated on a 1–5 semiquantitative scale in Figures S7 and S8 (^{footnote 4}). Although the lack of direct connection to spatial locations, geometry, or direct compositional information precludes much general insight, it is clear that the CJDS is extremely diverse on a textural level. A direct comparison of dike interior textures to elevation data and width (the only spatial data available in the WHT notebooks) does not result in significant correlations.

Dr. William H. Taubeneck documented ∼3% of the total CJDS dikes as containing some host-rock partial melt. We further investigated whether there was a relationship between any other variable in our database and thickness of the partially melted region. Figure S9 (^{footnote 4}) suggests a correlation between dike thickness and partial melt thickness. However, caution is required because we are not aware of what criteria WHT used to map partial melt. We do not attempt to compare with other studies of partial melt around the CRBG dikes (e.g., Petcovic, 2003).

Previous authors have identified distinct times at which joints and dikes interact (e.g., Delaney et al., 1986). This interaction can be interpreted as reflecting the ratio between the confining tectonic stress and the magmatic overpressure in concert with the preexisting fractures (e.g., Delaney et al., 1986; Jolly and Sanderson, 1995). Dike-joint interactions may also be superimposed on a larger overarching orientation (termed “zigzag” by Hoek, 1991). Dikes within the CJDS dike swarm appear to fit into this latter category. At the 1:24,000 map scale, dikes are linear features cutting across rock-type changes and even accreted terrane boundaries (Fig. 2A). However at the outcrop scale, some dikes appear to correlate in orientation with joints measured in granitoids (Fig. 11). This addresses the Taubeneck (1997) hypothesis, which posited that primary dike orientations were controlled by host-rock jointing. We instead interpret host-rock joints as influencing magma transport in dikes at an outcrop scale (less than tens of meters) only. Dike orientation at a map scale (1s–100s km) likely reflects regional tectonic or magmatic stresses at the time of emplacement (e.g., Corti et al., 2019).

The Chief Joseph dike swarm is the youngest and one of the best exposed magmatic plumbing systems of a large igneous province on Earth. This study synthesizes the unpublished dike mapping and observations made by Dr. William H. Taubeneck (Oregon State University; 1923–2016). Over the timespan of his career, WHT mapped 4279 individual dike segments with detailed notebook entries that document thickness, host-rock character, and qualitative textural observations. A database based on his work allows us to create the first holistic map of exposed dikes within the Chief Joseph dike swarm. We built on the WHT database by including dikes mapped by the Washington and Idaho Geological Surveys (Fig. 2A), supplemented by our own original observations and measurements. We also compiled the locations of all known mapped dikes related to the CRBG; those of the CJDS, Monument, and Steens swarms and the Ice Harbor dikes (Plate S1 [^{footnote 1}]).

The 1:24,000 scale map of dikes derived from the WHT database indicates that there are multiple and previously unrecognized systematic dike orientations within the CJDS (Fig. 2D). The dominance of a single NNW trend (Figs. 2D and 3A) suggests that throughout much of the active intrusive period, the least compressive horizontal stress was oriented WSW-ENE. Other orientation trends within the swarm largely postdate the main trend based on crosscutting relations.

On a sub-swarm scale, clusters of dike segments with similar orientations occur throughout the CJDS and point to a rich meso-scale structure of the swarm (Fig. 3). Clustered segments reinforce the dominance of the NNW trend in orientations, but overlapping clusters with different orientations also could indicate time-varying principal stresses. Clusters may reflect both dynamic dike-dike stress interactions as well as structural control of dike pathways by earlier structures. Further work is still needed to identify whether there is a clear progression in dike orientation.

On a dike-segment scale, dikes of the CJDS often segment and exhibit ∼1 m undulations in strike and width. When examined on the outcrop scale (1s–100s of m), dike segment orientations sometimes correlate with the orientation of joints within local granitoid rocks (Fig. 11). However, we believe joints are a second-order control on the form of dikes, behind regional-scale tectonic and magmatic stresses, reflected in the dominant NNW dike trend (Fig. 2D).

The number of dikes expressed regionally and the short time period involved have significant tectonic implications. We estimate a maximum observed dilation oriented E-W of 500 m with 150 m of dilation oriented N-S due to dike emplacement. If extrapolated regionally over areas where dikes are not well exposed, this increases to 2.7 km E-W and 1.3 km N-S dilation from CJDS. This dilation exceeds geologically observed amounts of strain in the proximal and largely concomitant northern Nevada rift (John et al., 2000). CJDS intrusion thus seems capable of transiently overcoming regional tectonic extension, consistent with time-varying variations in dike orientation established through dike crosscutting relationships (Fig. 9). Additionally, we calculate the upper bound volumes for the CJDS, ranging from 2.5 × 10²–6 × 10⁴ km³. This volume may contribute to part of the petrologically derived ∼160,000 km³ of intrusions associated with the erupted portion of the GRB.

We find a possible swarm-scale narrowing of dike segments with increased elevation and hence shallower paleodepth. This pattern is expressed over 1.5 km of paleodepth from the CRBG-granite contact exposed at ∼3000 m above sea level (Fig. 10). This trend could arise from viscous losses in magmatic overpressure with increasing distance from a deep magmatic reservoir; dike tipward thinning, or lateral dike propagation at a level of neutral buoyancy greater than ∼2 km depths.

The immensity of the CJDS dikes presents an intriguing challenge for CRBG eruption mechanics. The total volume of magma remaining in dikes implied by the observed CJDS is probably small compared to the erupted volume but is large enough and spatially distributed enough to affect crustal evolution and volatile budgets for the CRBG. Flow-rate calculations combined with the rarity of host-rock melting observations in the WHT data set suggest that the vast majority (perhaps as high as 97%) of exposed dikes were not long-lived transport structures and thus likely not feeders of surface flows.

A significant conclusion arises from the brief geochemical comparison of WHT dike samples with CRBG lavas: the dike swarm includes rocks that lie in non- or thinly populated regions of lava compositional space and extend the range of known CRBG chemistry (Fig. 7). This observation is consistent with only a subset of the mapped dikes feeding surface flows.

We view establishing the connection of the CJDS upwards to surface eruptions and downwards to magma storage zones and evaluating the mechanical consequences of CJDS emplacement into the crust at all scales as promising directions for future work. The quiescence of the rocks today hides the frenetic pace at which magma intruded this region, which is unique amongst dike swarms on Earth in its scale of exposure and documentation, quality of exposure at a range of scales, and rich history of research. Apart from being an exquisite testament to a tenacious field geologist, we hope that the legacy data set of William H. Taubeneck, in turn, provides a baseline for future research.

We would like to first thank William H. Taubeneck for his years of diligent mapping and work on the geology of northeastern Oregon and study of the Chief Joseph dike swarm. We also appreciate the aid provided by the Whitman College and University of Oregon libraries in scanning the Taubeneck collections. Our thanks also go to the three years of University of Oregon Field Camp students who helped us collect important primary observations. Dr. Karlstrom acknowledges funding from National Science Foundation awards EAR-1547594-001 and EAR-1848554. This manuscript was much improved by the feedback and conversations with Drs. Eugene Humphreys and Ray Weldon. Connor Mullady and Scott Boroughs re-analyzed the old WHT samples in the Washington State University GeoAnalytical Lab. We appreciate helpful and insightful reviews from Drs. Vic Camp and Francesco Mazzarini, Science Editor Shanaka de Silva, Associate Editor Valerio Acocella, and one anonymous reviewer.