Large‐Volume and Shallow Magma Intrusions in the Blackfoot Reservoir Volcanic Field (Idaho, USA)

The Blackfoot Reservoir volcanic field (BRVF), Idaho, USA, is a bimodal volcanic field that has hosted silicic eruptions during at least two episodes, as recently as 58 ka. Using newly collected ground and boat‐based gravity data, two large negative anomalies ( −16 mGal) are modeled as shallow ( <1   km) intrusions beneath a NE‐trending alignment of BRVF rhyolite domes and tuff rings. Given the trade‐off between density contrast and model volume, best‐fit gravity inversion models yield a total intrusion volume of 50−120 km3 ; a density contrast of −400 kg   m−3 results in two intrusions, each ∼9 km ×4.5 km and about 0.5 km thick, with cumulative volume of 100   km3 . A network of 340°−360° trending faults lies directly above and on the margins of the mapped gravity anomalies. Most of these faults have 5−10 m throw; one has throw up to ∼50 m. We suggest that the emplacement of shallow sill‐like intrusions produced this fault zone and also created a ENE‐trending fault set, indicating widespread ground deformation during intrusion emplacement. The intrusions and silicic domes are located 3−5 km E of a regional, 20 mGal step in gravity. We interpret this step in gravity as thickening of the Upper Precambrian to lowermost Cambrian quartzites in the Meade thrust sheet, part of the Idaho‐Wyoming Thrust Belt. Silicic volcanism in the BRVF is a classic example of volcanotectonic interaction, influenced by regional structure and creating widespread deformation. We suggest volcanic hazard assessments should consider the possibility of large‐volume silicic eruptions in the future.

• Large-amplitude gravity anomalies are mapped in a combined ground and boat-based gravity survey in the Blackfoot Reservoir volcanic field, Idaho (BRVF), adjacent to young (1.5 Ma, 58 ka) topaz rhyolite domes and tuff rings within a Quaternary basaltic volcanic field • Best-fit 3D inversion of the gravity data, constrained by density contrast estimates and excess mass calculations, indicates the presence of two bodies with thick sill-like shapes in the uppermost crust, with cumulative volume of ∼100 km 3 and volume uncertainty in the range The Blackfoot Reservoir volcanic field (BRVF), located in the northeast Basin and Range of the western USA (Figure 1), is a bimodal volcanic field (McCurry & Welhan, 2012). We use newly collected groundbased and boat-based gravity data to investigate the potential for shallow intrusions associated with an alignment of five silicic domes and explosion craters, erupted approximately 58 ka, in an area called the Central Dome Field (CDF) located within the BRVF (Figure 2a). The CDF includes a network of N to NNW-trending surface faults that are unique to the region in their variable along-strike displacement and en echelon, corrugated map pattern (McCurry & Welhan, 2012;Polun, 2011). These features suggest that these are young normal faults (Ferrill et al., 1999), similar to those produced by volcanotectonic interaction mapped in other volcanic fields Bursik & Sieh, 1989;Garibaldi et al., 2020;Gottsmann et al., 2009;Mazzarini et al., 2004;Tuffen & Dingwell, 2005) (Figures 2b and 3), inspiring us to further evaluate the potential for shallow intrusions.
We use 3D gravity models to explore the potential subsurface geometries that create the observed gravity anomalies we map in the CDF. The models are calibrated with the density of nearby silicic domes and with an excess mass calculation. We estimate the volumes of the inferred intrusions and the domes to constrain the intrusive to extrusive volume ratio. The locations and displacements of faults (McCurry & Welhan, 2012;Polun, 2011) are compared with the boundaries of the modeled intrusions (Figures 2a,2b,and 3). Our results suggest that potential future silicic activity may involve comparatively large volumes of silicic magma and could be accompanied by widespread surface deformation. Results also suggest that regional tectonic structures may influence magma ascent and accumulation in the shallow crust, as found in other volcanic systems (Acocella & Funiciello, 1999;Bacon et al., 1980;Deng et al., 2017;White et al., 2015).

Overview of BRVF Geology
The BRVF lies in the transition between the Intermontane Seismic Belt and a seismically quiescent region that includes the Eastern Snake River Plain (ESRP) (Anders et al., 1989). The BRVF sits roughly 200 km from Yellowstone and has experienced younger eruptions than the most recent eruptive lava flows at Yellowstone (Eaton et al., 1975). This distributed volcanic field comprises Quaternary scoria cones, basalt flows, rhyolitic domes, and tuff rings ( Figure 3). There are three rhyolitic domes at the southern end of the Blackfoot Reservoir, named China Hat, China Cap, and North Cone. These three domes and nearby tuff rings make up a NE-trending volcano vent alignment that defines the CDF (Figure 2b). The base of the China Hat and China Cap domes are primarily block and ash flows with surge deposits exposed in a quarry at the base of China Hat dome. The craters of two tuff rings, Burchett Lake and Gronewell Lake, are filled with water. These tuff rings have low outer slopes typical of surge deposits associated with phreatomagmatic eruptions (Figure 2b). The China Cap dome has been dated using 40 Ar/ 39 Ar, yielding an age of 58 ka (Heumann, 2004). Sheep Island lies on the western side of the Blackfoot Reservoir and is dated to part of a prior eruptive episode roughly 1.5 Ma and provides evidence for silicic volcanism in the BRVF to be recurring rather than a singular event.
The basaltic lavas of the BRVF erupted from low scoria cones and fissures. Basalt lava flows reach a thickness of 290 m in the CDF, where they surround the silicic vents and cap the underlying geology as a continuous lava flow field. Basalt eruptions in the BRVF have poor age constraints. Some of the lavas from the BRVF flowed out to the southwest into Gem Valley (Figure 1). These have been dated radiometrically between 100 and 25 ka (McCurry et al., 2011). Basalt vent alignments also occur in Gem Valley.
Mapping of the surrounding bedrock geology reveals several generations of faults including NW-trending, SW-dipping thrust faults of the Idaho-Wyoming Thrust Belt (Figures 2 and 3) formed during the Jura-Cretaceous Sevier Orogeny (Armstrong & Oriel, 1965;Dixon, 1982). NW-trending normal faults, perhaps representing two phases of late Tertiary extension, overprint these older faults. In addition to these older structures, there is a third set of distinctive normal faults (Polun, 2011) (Figures 2 and 3) that are only found within the BRVF. We evaluate the origin of these latter faults and their relationships to silicic volcanic vents in light of gravity anomalies and models, described below.  Mabey and Oriel (1970) first identified negative gravity anomalies in the CDF, which they interpreted as shallow sedimentary basins. Using the same gravity data set, Leeman and Gettings (1977) interpreted the gravity anomalies in the CDF to be related to a large silicic magma body (∼330 km 3 ) in the upper 6 km of the crust. Their model of the gravity anomalies in the CDF is consistent with the spatial association of the anomalies with young silicic domes. Also, this interpretation is consistent with prominent gravity anomalies associated with shallow intrusions elsewhere (Battaglia et al., 2003;Blakely, 1994;Bott & Smithson, 1967;Finn & Williams, 1982;George et al., 2016;Miller et al., 2017;Paulatto et al., 2019). We provide newly collected ground-based and boat-based gravity data that further constrain these anomalies and their relationship to Quaternary volcanoes and faults.

Gravity Data Collection and Processing
New gravity data were collected broadly throughout the BRVF, with higher density sampling in and around the CDF to identify the shapes of the anomalies. These data were merged with the regional database (Keller et al., 2006), consisting almost entirely of data collected by the USGS, including survey data collected by Mabey and Oriel (1970). In addition to ground-based data, we collected boat-based gravity data over the A total of 460 new ground-based gravity measurements were made with a Burris gravimeter (B-38) with measurement precision of approximately 0.003 mGal. Station location was determined using a Trimble R10 and CenterPoint RTX service, which has a horizontal precision of 3 − 5 cm and a vertical precision of 7 − 10 cm (Glocker et al., 2012). After correcting for an instrument drift of ± 0.025 mGal/day, the uncertainty on our gravity measurements is ± 0.03 mGal. The gravity base station is in the town of Soda Springs and the same base station was used throughout all of the campaigns. This allowed us to make multiple base readings each day of the survey to accurately capture the instrument drift, which is quite linear for this instrument.
Ground-based gravity data reduction included tidal, latitude, atmospheric mass, free-air, spherical cap Bouguer and terrain corrections (White et al., 2015). These corrections were applied to the new data and to the drift-corrected regional data from the USGS to achieve consistency among gravity data from different sources. The terrain correction was applied in two parts. An inner correction used a 10 m DEM with 20 km radius about each gravity station, and an outer correction used a 30 m DEM with 167 km radius about each station. The DEM data used for the terrain corrections were obtained from the USGS National Elevation Database (NED), and a density of 2,670 kg m −3 was used for Bouguer and terrain corrections as it is generally accepted as the average density of crustal rocks (Hinze, 2003). Gravity was remeasured at several USGS gravity station locations to use as tie-in points, similar to the procedure in Deng et al. (2017).  (Armstrong & Oriel, 1965). The Hubbard 25-1 borehole is represented by the green star. Red triangles show basaltic vents.
The ground-based gravity data reveal large amplitude (∼21 mGal) negative anomalies in the CDF with a gravity gradient under the reservoir (Figure 4). We collected over 14,000 data points with a Dynamic Gravity Systems (DGS) Marine Gravity Sensor (AT1M) on a pontoon boat to define the shape and gradient of the gravity anomaly in the reservoir (Figure 5b). This gravimeter is gimbaled to compensate for the accelerations imposed by the motion of the boat. The same corrections made to the ground-based data were applied to the boat-based data, with additional corrections accounting for the motion of the gravimeter. The Eötvös correction was applied to account for the velocity of the boat as it adds or subtracts to the tangential velocity of the gravimeter relative to the rotational axis of Earth, and the acceleration of the platform the gravimeter rests on was accounted for in the inertial reference frame of the vessel (Telford et al., 1990). A correction was made for the mass of water in the reservoir, although this is found to have trivial impact as the reservoir is <10 m deep and changes depth very gradually (Wood et al., 2011). The velocity and acceleration of the vessel were obtained through the differentiation and double differentiation of the GPS position, respectively.
The boat-based data were sampled at a rate of 1 Hz on a continuously moving platform, leading to a higher spatial density of measurements on the reservoir compared to the ground-based measurements. Including  Oriel and Platt (1980), showing that the Quaternary basalts cover the valley floor and flowed toward the town of Soda Springs to the south and Gem Valley to the southwest. The faults in the BRVF show a distinctly different trend/orientation relative to the bedrock faults in the adjacent ranges.
all of the boat-based data in our gravity model would cause the region beneath the reservoir to be over-constrained leaving the more sparsely sampled ground-based regions to be comparatively under-constrained and less significant in the gravity model. Consequently, the boat-based data were sampled every 100 m along the survey track lines to mitigate over-constraining the region beneath the Blackfoot Reservoir during the inversion.
The combined ground-based and boat-based data were further filtered to include only a 780 km 2 area (3,126 measurements), centered on the two negative CDF gravity anomalies (Figure 5b). This filtering helps to identify longer wavelength, regional signals that underlie the negative anomalies in the BRVF and to separate these shorter wavelength gravity anomalies from the regional gravity, as described in the next section. Both the entire data set and the grid of sub-sampled data used to model the anomalies are provided in the supplementary material.
Overall, the data coverage defines the two gravity anomalies in the CDF, including the continuation of one of the anomalies beneath the reservoir as discovered and verified by the boat survey. Gravity stations are most numerous over these anomalies and in the high-gradient areas adjacent to the anomalies. Therefore,   (Evans, 2010), with faults, domes, and vents. Normal faults are marked by black lines with throw markers; ENE trending faults southeast of the rhyolitic domes (red patches) are black lines without throw markers. Basaltic vents are red triangles. The Hubbard 25-1 borehole (green star, Figure 6) is located just south of China Hat dome. The map region is constrained to the data bounds used for the inversions. (a) Reference map centered on the BRVF used as the bounded region for the inversion, (b) complete Bouguer gravity with all gravity stations represented by colored circles, (c) regional and (d) local gravity anomalies with the gravity stations used in their inversions respectively. All the gravity stations have a 1 km radius mask to highlight the best constrained areas. Blue grid lines show the prisms boundaries used in the respective inversions. Prisms for the regional model (c) are 4 × 4 km and extend slightly past the data bounds to minimize edge effects; prisms for the local model (d) are 2 × 2 km. the gravity anomalies are well resolved by the gravity station distribution. One exception is the SE corner of the southern negative anomaly. In this area we must rely on the older USGS gravity data set, as access was not permitted to these lands by the current property owner. Nevertheless, dense sampling up to the edge of the property and the older USGS data reasonably constrain the gravity gradient, just not with the resolution obtained elsewhere on the map.

Isolation of the CDF Gravity Anomalies
Gravity anomalies arise from a combination of broader regional effects of the basement structure and shorter wavelength anomalies produced by local mass variations in the shallower subsurface. Separating the local gravity anomalies from the regional gravity signal is paramount to interpreting and modeling the gravity data. The complete Bouguer gravity map of the CDF (Figure 5b) includes two distinct, negative gravity anomalies with magnitude of approximately −21 mGal. These short wavelength anomalies lie within a regional gravity anomaly, with high amplitude positive values (20 mGal) to the west and low amplitude negative (−5 mGal) values to the east ( Figure 4). The regional variation does not correlate with the topography, and the transition between the positive and negative values happens over a relatively short distance (∼8 km). This gradient is not linear, but shows a step in the regional gravity that is located 2 − 3 km west of the rhyolite domes in the CDF ( Figure 5b).
The regional gravity trend was isolated by removing data more negative than a threshold value of −6 mGal, which was chosen by graphical separation of the local minima within the regional trend ( Figure 5c). The filtered data that were removed are interpreted to be the local gravity anomalies. The threshold value used to separate the regional anomaly from the local is subtracted from the local data and these data are contoured ( Figure 5d). The filtered local gravity anomaly has an amplitude of approximately −15 mGal, with clear separation from other sources of anomalous gravity. Adding the two maps (Figures 5c and 5d) gives the original gravity map (Figure 5b).
The regional, long-wavelength gravity anomaly ( Figure 5c) shows a large amplitude positive anomaly (20 mGal) over the range between Gem Valley and the BRVF. A cross-sectional profile from Dixon (1982) (his number 17) depicts the west-dipping Meade thrust fault cutting and displacing the contact between the Precambrian and Cambrian (1−3 km depth). This displacement shallows and thickens quartzites beneath the range on the western edge of the BRVF. We suggest that the observed regional gravity step correlates to the approximate eastern limit of the quartzites that are displaced in the Meade thrust fault.
The local gravity anomalies have elliptical shapes, each striking NW− SE. The two negative anomalies are separated by a saddle of higher gravity values ( Figure 5d). The domes and tuff rings lie within and near this saddle. The volcano vent alignment is nearly orthogonal in trend to the long-axes of the negative anomalies. The faults in the BRVF appear to wrap around the negative anomalies on the west side of China Hat dome and the western margin of Blackfoot Reservoir (Figure 5d).

Constraints on the Local Gravity Model
The two negative CDF gravity anomalies ( Figure 5d) represent a mass deficit. We calculate the mass deficit, Δ , using Green's function (Parker, 1974): where Δ ( ) is the gravity anomaly, and are the number of grid points in the (easting) and (northing) directions, respectively, and Δ and Δ is the grid spacing (500 m) in the and directions. This integration of the detrended gravity data gives a mass deficit of −3.5 × 10 13 kg. For a reasonable range of density contrasts, the mass deficit calculation shows that the causative body of these anomalies is of order of one hundred cubic kilometers of material.
Hand samples of rhyolite from the China Cap dome yield unsaturated bulk rock densities of 1600 − 1800 kg m −3 . The Nettleton and Parasnis approaches to modeling bulk density from gravity profile data (Agustsdottir et al., 2011;Nettleton, 1939;Parasnis, 1952;Saballos et al., 2013) yield a bulk dome density of about 1,700 kg m −3 for China Cap dome, which is consistent with bulk silicic dome densities determined using the same methods elsewhere (Agustsdottir et al., 2011). We assume that the density contrast between intrusive silicic rocks and the crust is not as large as the density contrast between the rhyolite dome and the crust, but it may approach this value. Additionally, density estimates of A-type granophyres and rhyolite intrusions are as high as 2400 kg m −3 (Lowenstern et al., 1997). The Hubbard 25-1 Borehole (Figure 2b), drilled in 1983, provides constraints on the density and lithology of the country rock within the upper crust of the BRVF (Polun, 2011). The well is located approximately 1.5 km south of China Hat and approximately 1 km west of the edge of the southern negative gravity anomaly ( Figure 5b). The compensated neutron lithodensity logs contain data that constrains the bulk density as a function of depth within the borehole. The range of densities within the log spans from 2600 − 2800 kg m −3 with an average density over the entire 2 km section of 2,700 kg m −3 ( Figure 6). The lithology within this well alternates between basalts, siltstones, and shales near the surface to interbedded limestones, sandstones, and shales at depth. The thickness of basalts in the uppermost part of the log is approximately 290 m including scoria layers, constraining the thickness of BRVF basalts. We were unable to determine from the logs if the deeper basalts (750 and 1,100 m) are extrusive or intrusive. Nevertheless, we are confident that igneous rocks are present at these depth intervals.
Given a mass deficit of −3.5 × 10 13 kg, for density contrasts −800 to −300 kg m −3 , the causative body has a volume range of 50 − 120 km 3 . This range of density contrasts is used in our gravity inversion models and our model results are compared with this range of volume estimates.

Gravity Modeling of Regional and Local Anomalies
Inverse modeling is used to deduce subsurface structure both for regional and local anomalies (Figures 5c  and 5d). Our modeling approach first discretizes the subsurface into a grid of vertical-sided rectangular prisms (i.e., the blue grids in Figures 5c and 5d). We assume a constant density contrast between all prisms and the surrounding bedrock, but the magnitude of this density contrast is solved during inverse modeling of the gravity data.

Inversion Procedure
Two inversion procedures are used, one to model the regional signal and one for the local anomalies. Regional inversion modeling assumes a single bottom depth for all prisms, while local inversion modeling uses unique top and bottom depths for each prism. Inputs to the inversion include a range for each adjustable parameter value (depth-to-bottom, depth-to-top, density contrast). Both inversions initialize multiple sets of initial parameter guesses, drawn from input ranges specified in a configuration file. The total number of parameter sets is one more than the total number of modifiable parameters. The local inversion model has 391 independent model parameters, resulting in the initialization of 392 unique sets of randomized parameters; the regional inversion model has 58 independent model parameters, resulting in the initialization of 59 unique sets of randomized parameters.
The inversion process adjusts and tests these parameter combinations, using a calculated solution for the gravity due to a prism. The gbox solution for gravity (Blakely, 1996), written in C for speed, is used as the forward model. The gravity anomaly associated with each prism is summed across the map area and then compared with observed gravity values interpolated on to a grid. Interpolated and gridded gravity values are used because of variability in the density of gravity measurements across the region and to speed calculations. The grid size for the inversion process is selected by experimentation to minimize the number of model parameters and to best resolve the subsurface structure. Modeling a large number of small prisms often results in an awkward prism solution that requires additional smoothing, which does not necessarily improve the model (White et al., 2015). Our modeling attempts using a large number of small prisms created unrealistic bumps and rapid changes in prism thickness, resulting in an unrealistic model geometry given the relatively smooth variation in the observed gravity.
The downhill-simplex optimization algorithm (Nelder & Mead, 1965;Press et al., 2007) is used to resolve and identify a best set of model parameters based on a goodness-of-fit test designed to minimize the residual error between the measured data and the calculated solution. We use the root-mean-squared error (RMSE) for this goodness-of-fit test. Typically, 100, 000 − 200, 000 forward solutions are calculated to find a best-fit model. Multiple simulations are completed by varying the random seed and prism boundaries to fully explore the model parameter space and to identify local minima.  Figure 2). The average host rock density through the upper 2.5 km in the BRVF is 2,700 kg m −3 , a higher than average density that adds to the density contrast causing the negative CDF gravity anomalies. The unit from ∼500− 600 m was characterized as anhydrite in the log report, although the low density of this unit conflicts with the normally high density of anhydrite (Robertson et al., 1958). Other units are high density sedimentary rocks consistent with the Paleozoic section and some basalt sequences that may be intrusions within this sequence, capped by Quaternary basalt.

Regional Model
The model of the regional gravity field (Figure 5c) is based on the interpretation that a thickening of Precambrian quartizites in the Meade thrust fault exists near the western edge of the BRVF (Dixon, 1982). The prism size used for the regional model is 4 × 4 km, due to the more widely-spaced gravity data to the west of the BRVF. We model the regional data with a flat-bottomed geometry to more closely emulate the thickening of quartzites on the west side of the BRVF. The modeled density contrast ranges from 0 to 150 kg m −3 and the modeled depth range for the quartzite contact is 0.5 − 12 km. The model prisms extend slightly beyond the data boundaries to resolve edge effects and better constrain the gravity anomalies at the edges of the model area (Figure 5c). Figure 7 shows the geometry of the best-fit inversion model for the regional gravity data. The depth-tobottom is 8.1 km; all models solved for a density contrast around 150 kg m −3 . The average depth-to-top on the western margin of the region is ∼2 km, which is in agreement with the range from Dixon (1982) for the Inversion results for the regional gravity anomaly. The top perspective image depicts the CDF over the extent of the prisms for the inversion of the regional anomaly. The centers of the prisms are represented by circles that are colored and contoured by the depth to the tops of the prisms. The bottom depth of this model is uniform at 8.1 km and the model density contrast is 150 kg m −3 . The bottom plot is a 3D perspective mesh of the tops of the prisms and is colored by depth-to-top. This model shows that a thickening of high density quartzites associated with thrust faulting is a possible cause of the regional anomaly.
depth to the Precambrian-Cambrian contact (between 1.5 and 3 km). The regional model shows that the quartzites are thickened by 6 km, on average, near the range on the western edge of the BRVF, and that the Precambrian-Cambrian contact sits at roughly ∼8 km depth in the area of the local anomalies of the CDF. The shallowest prisms in the model are in the southwestern region of the model where it reaches a depth of ∼ 650 m where the highest gravity values are located (∼20 mGal). The regional model is not able to reproduce the highest gravity values (>18 mGal) without increasing the density contrast, but a higher density contrast does not agree with known densities of quartzite. The model suggests that the regional step in the gravity field is related to the approximate eastern limit of the thickening quartzites in the Meade thrust sheet, but the story is likely more complex.

Local Model of the Igneous Intrusions
Inversion models of the local CDF gravity anomalies (Figure 5d) are constructed using a wide range of potential density contrasts (−100 kg m −3 to −900 kg m −3 ). The minimum value for the depth-to-top parameter is 250 m, based on the approximate thickness of the basalt section (McCurry & Welhan, 2012). This lithologic and mechanical contrast is assumed to introduce a mechanical and compositional boundary that would limit the depth to the top of the intrusions (Kavanagh et al., 2006;Richardson et al., 2015;Wetmore et al., 2009). The maximum value for the depth-to-bottom parameter is constrained to 3 km. Maximum prism depths deeper than 3 km tend to produce anomalies of longer wavelength than the observed anomaly.
All best-fit models show two compact bodies in the shallow (<1 km) subsurface that thin toward their margins, giving them a sill-like geometry (Roman-Berdiel et al., 1995); the two bodies have thin or absent prisms between them. Density contrasts between −800 and −500 kg m −3 tend to produce geometries with more variation in depth to top of the bodies, a laccolith shape, while density contrasts between − 300 and −500 kg m −3 tend to produce geometries with more variation in depth to bottom and bodies with flatter tops, a lapolith shape. As in all gravity models, there is parameter compensation in the tradeoff between density contrast and volume. For example, increasing the density contrast can result in thinner prisms on average, and conversely, decreasing the density contrast can result in thicker prisms. We tested and compiled best-fit models by imposing limits on the density contrast to evaluate the tradeoff between volume and density contrast of the model space. Some of these model results did not have low RMSE. Larger density contrast results in a deeper average depth of the body, but all are relatively shallow (average depth ≤ 1 km). Figure 9 shows the solutions for 17 simulations, each testing 100, 000 − 200, 000 parameter combinations. This plot illustrates the tradeoff between density contrast and volume (Blakely, 1994). Solutions have density contrasts between −800 and −400 kg m −3 and agree with: (a) lithology observed in the Hubbard 25-1 borehole, (b) dome density determined from China Cap hand samples and Parasnis/Nettleton density analyses (Nettleton, 1939;Parasnis, 1952), and (c) volume estimates from mass deficit. A range of reasonable solutions with nearly identical RMSE occur between density contrasts of −600 to −400 kg m −3 . These solutions give a range of volume estimates from ∼60 to ∼100 km 3 . The minimum volume of the anomalous mass is ∼50 km 3 with a maximum density contrast of approximately −800 kg m −3 . The maximum volume of ∼120 km 3 is obtained with a density contrast approximately −300 kg m −3 , acknowledging that the RMSE is higher for this low density contrast model. In all models the northern body is larger than the southern body. For example, at −400 kg m −3 density contrast the volume of the northern anomaly is approximately 60 km 3 and the volume of the southern anomaly is approximately 40 km 3 .

Modeling the Gravity Anomalies as Shallow Intrusions
The new gravity data, combined with previous surveys, identifies two large negative anomalies. The addition of boat-based gravity data constrains the western margin of the northern gravity anomaly, which resides largely under the Blackfoot Reservoir. Based on these data and models, we suggest that the large negative gravity anomalies within the CDF are due to high-level silicic intrusions rather than due to a sedimentary basin, as inferred by Mabey and Oriel (1970). If the anomalies were produced by sediments, the basin would be thickest toward the center and the anomaly would have low gravity gradient near its center (Gimenez et al., 2009). Instead, the anomalies show short-wavelength variation where they have the largest negative values. These short-wavelength anomalies indicate that the causative body is actually closer to the surface near the centers of the gravity anomalies. We tested the sedimentary basin model and found poor fits (high RMSE) to the observed gravity data, especially in the center regions of the isolated negative gravity anomalies where the amplitude of the anomalies is high. It is particularly difficult to model basin geometries that create a narrow divide between the two isolated depocenters.
Geologic data support the interpretation that the gravity anomalies are related to igneous intrusions rather than to sedimentary basins. One key observation is from the Hubbard 25-1 exploration log (Polun, 2011). The presence of anhydrites in the upper 700 m suggests that the area of the CDF was submerged and gradually infilled by sediments eroded from the adjacent ranges. However, this section is thin (∼400 m) and has a small density contrast indicating that it is unlikely the negative gravity anomalies are related to this stratigraphic sequence. Additionally, we note the anhydrite unit in the well log ( Figure 6) is logged as a lower density unit, which is inconsistent with the high density of anhydrite (Robertson et al., 1958). It is possible the anyhdrite unit is actually misidentified in the log. The rest of the section is dominated by a passive margin sequence characteristic of the Paleozoic section. Figure 8. Example inversion results for the local gravity anomalies. The modeled density contrast is −400 kg m −3 ; the deepest prism extends to a depth of 2.9 km. Thickness contours of the modeled prism geometry (a) are plotted over a 10-m hillshade DEM with faults, vents, and domes superimposed. Model prisms with thickness 100 m, are outlined with blue squares that underlay the thickness contours. A 3D perspective of the prism geometry with 2.5 times vertical exaggeration (b) illustrates the separation between the two distinct bodies modeled by the inversion. Basaltic vents and rhyolitic domes are represented by red and black triangles respectively; faults are marked by black lines with fault throws; location of the Hubbard 25-1 borehole, detailed in Figure 6, is depicted by a green star (a) and green cylinder (b). Animations of the 3D rendering can be found in the supplementary material.
There is an absence of clear basin-bounding normal faults on the eastern and western margins of the BRVF, whereas sedimentary basins in the region have clear basin-bounding faults. The west margin of the modeled intrusion coincides with a west dipping fault with the largest vertical offset (50 m) observed in the BRVF. This sense of offset is concurrent with deformation during the emplacement of shallow intrusions (Acocella, 2000;Acocella et al., 2002;Castro et al., 2016). We note that the sense of offset is opposite of that which would be expected if the fault bounded a sedimentary basin. There are plenty of basins in the region, Gem Valley for example, but all are elongate parallel to basin-bounding faults and none of them exhibit this pattern of faulting.
Density contrast between the interpreted intrusions and the country rock is a source of uncertainty. The density of the country rock is constrained to be approximately 2,700 kg m −3 . The densities of dome rocks, 1,700 kg m −3 , are likely too low and produce too high a density contrast compared to high-level intrusive equivalents to these dome rocks. Rhyolite melt densities are typically 2,350-2,400 kg m −3 (Bachmann & Bergantz, 2004), which would produce a density contrast of approximately −300 to −350 kg m −3 . Granites are created from rhyolite magmas in the midcrust through crystallization of dense mineral phases, filter pressing and compaction, all of which leaves a lower density residual melt that can ascend to high crustal levels or erupt (Bachmann & Bergantz, 2004). These bodies, interpreted to be high-level intrusions, are also well below saturation pressures for volatiles in silicic magmas, and so may be porous and may be fractured Figure 9. The trade-off between density contrast and volume of bodies associated with the local gravity anomalies (Figure 8) is illustrated using 17 different inversions. Each circle represents an inversion result; the size/color of the circle corresponds to the goodness-of-fit (RMSE) of the inversion. Inversion results give a minimum intrusion volume of 50 km 3 with a maximum density contrast of −800 kg m −3 . A range of reasonable solutions between −600 and −400 kg m −3 that have respective volumes between approximately 60 and 120 km 3 is identified by the blue box. and altered during cooling. Both of these processes result in lower bulk rock density. For example, 10% saturated bulk porosity in a rock of nonporous density of 2,350 kg m −3 yields a bulk density of 2,260 kg m −3 , or a density contrast of −440 kg m −3 . Density contrasts of around −600 kg m −3 are used to model gravity anomalies associated with other high-level intrusions (Acocella, 2000;Miller et al., 2017).
We suggest a reasonable range of density contrasts between the intrusions and the country rock is −600 to −400 kg m −3 . All models in this range of density contrasts produce two elliptical shaped bodies each with approximate map dimensions of ∼ 9 km × 4.5 km. Altering the density contrast in this range results in a thickening or thinning of the intrusions while the horizontal dimensions remain relatively constant (Figure 10). This range of density contrasts corresponds to cumulative intrusion volumes of 60 − 100 km 3 . A density contrast of −400 kg m −3 yields an intrusion volume of approximately 100 km 3 .
Both gravity anomalies, and by inference the intrusions, are slightly elongate NW, perpendicular to the NE (approximately 35 • ) alignment of silicic domes (Figure 5d). This geometry is consistent with the high-level intrusion model proposed by Vigneresse et al. (1999). In the absence of substantial volume of intrusion, the unperturbed stress state in the region is extensional, with 1 vertical and equal to lithostatic pressure in magnitude. A fracture or dike will propagate vertically and perpendicular to the least principle compressive stress, 3 . From the vent alignment we infer that 3 is oriented approximately 125 • . As the intrusion shallows, the magma pressure exceeds the lithostatic pressure causing a stress rotation, with 3 becoming vertical, resulting in horizontal intrusion. 2 becomes oriented approximately 125 • and 1 approximately 35 • , allowing the intrusion to grow faster in a NW-SE direction, perpendicular to the trend of the vent alignment.

Emplacement Related Deformation
The coincidence of the edges of the negative gravity anomaly with dramatic, if relatively small displacement faults points to volcanotectonic interaction during intrusion and silicic dome eruptions (Bursik & Sieh, 1989;Bursik et al., 2003). The faults in the BRVF extend from just north of the town of Soda Springs through the Blackfoot Reservoir, only cutting through bedrock at the surface near the southern end of Pelican Ridge (Figure 2a). While Polun (2011) placed the eastern limit of the rift zone at the discontinuous Hole in the Rock-China Hat fault, we believe, based on topographic data available through the Idaho LiDAR Consortium (Figure 10), that the eastern margin of the rift is an unnamed fault located along the western slopes of the Fox Hills extending north to the east of the Blackfoot Reservoir (Figure 2). The maximum E-W width of the faulting in the BRVF, at the latitude of China Hat, is ∼10.7 km. The faults in the BRVF are primarily NNW to NNE-trending and exhibit both east and west dips.
The western portion of the fault system in the BRVF includes a prominent nested graben trending N to NNW with the most topographically well-defined portion located just west of the rhyolite domes ( Figure 10). The graben is bounded on the west by the east-dipping Government Road Fault, which has a prominent scarp that is as much as 50 m high. The Government Road Fault is flanked on its west in its central portion by two additional east-dipping faults with scarps as large as 15 m (Figures 2 and 10). The eastern side of the graben is defined by the west-dipping Hole in the Rock and China Hat faults, which appear to be separated by a small left step just north of the China Hat dome (Figures 2 and 10). The graben appears to be floored by a loess-covered surface that is composed of the lavas from several basaltic vents including Red Mountain. The surface steps down 100 m from west to east across a series of east and west-dipping faults creating narrow (∼50 −150 m) full and half grabens separated by relatively broad (∼250 − 750 m) horsts. Throughout the broader graben the surface is typically flat or dipping slightly (<3 • ) east, a slope that appears to have been, at least in part, present before the youngest phase of faulting based on profiles outside the graben to the north and south.
Polun (2011) estimated horizontal extension across the graben from fault displacement and dip. These estimates suggest that the portion of the horst and graben system most proximal to the CDF has the largest magnitude of horizontal extension ranging between 75 and 200 m, depending on the fault dips. The total extension is taken to be a minimum because the estimates did not include all of the faults on the eastern extent of the fault system. The estimates based on minimum extension (i.e., fault dip of 70 • ) indicate increases from single digits to >50 m over a distance of 4 − 5 km on either side of the CDF. Based on these data, it appears that extension in the BRVF is greatest adjacent to the gravity anomalies and silicic domes, consistent with faulting during emplacement and/or draining of the intrusions.
A set of ENE-trending faults are only found directly overlying the intrusions, especially SW of China Hat dome. These faults appear to be unrelated to the normal tectonic setting of the BRVF. Instead, these faults may have formed during uplift and possibly deflation associated with the intrusions, perhaps associated with the extrusion of magma at the nearby domes (Figure 5d). This ENE-trending fault set is far less pronounced than the other faults in the BRVF (Figure 2b). The average throw across faults in this set is 1 − 2 m with a maximum of ∼10 m. Most of the faults are north dipping with the exception of one in the northern third of the set and the three southern-most faults. Acocella and Funiciello (1999) show that roof lifting associated with the emplacement of a laccolith is viable in producing significant uplift over the intrusion as well as faulting at the margins of the intrusion. We suggest that the pattern of diffuse faulting at the surface is associated with the emplacement of the modeled intrusions and draining of the shallow magmatic system during eruption of the CDF rhyolite domes. The highly faulted graben on the west end of the CDF has the greatest extension and lies on the margin of the modeled intrusion geometry. This shows a spatial correlation with the margins of the intrusion and the greatest structurally accommodated extension (Spinks et al., 2005). The amount of horizontal extension that is accommodated is at minimum 75 − 200 m in the CDF. Castro et al. (2016) has shown that shallow (20 − 200 m), rapid intrusion of laccoliths can produce large uplift (>200 m) and deformation at the margins of intrusion. In the BRVF, we observe the highest magnitude of faulting near the CDF and gravity anomalies with waning surface deformation north and south of the Figure 11. Synthesis of data and the interpretation of our model. (a) The map shows a perspective of the CDF obliquely perpendicular to the strike of the gravity anomalies, parallel to basaltic vent distribution, and is overlain by the complete Bouguer gravity anomaly as well as the regional and local faults, vents, and domes. The cross section line AA′ runs SW− NE from Gem Valley to the Meade Thrust. (b) The complete Bouguer gravity anomaly and elevation along the profile line shows the gravity high on the western side of the BRVF and the two gravity lows in the CDF that are separated by a saddle, or relative gravity high bounded by two gravity lows. The cross section (c) illustrates the schematic interpretation of the gravity anomalies in the BRVF and structural geology of Dixon (1982). The gravity high on the west side of the BRVF correlates to the relative shallowing of Precambrian/Proterozoic quartzites via the Meade Thrust (Figure 7), and the gravity lows in the CDF are interpreted as shallow rhyolitic intrusions associated with the domes at the surface, and faulted, uplifted topography ( Figure 8). The source of the shallow intrusions is inferred to be much deeper (∼ 14 km) and likely related to the source of the rhyolitic magmas (dashed lines) (McCurry et al., 2015). gravity anomalies. Our model suggests that shallow silicic intrusions were emplaced, uplifted the BRVF and generated ancillary networks of faults similar to the Cordón Caulle (Castro et al., 2016).
Overall, the CDF and its twin gravity anomalies are closely associated with faulting on several scales. The area of the CDF is marked by two negative gravity anomalies, interpreted to be sill-like intrusions, and by faulted topography (Figures 11a-11c). These faults wrap around the two gravity anomalies, especially on the west side of the reservoir, a fault pattern that is consistent with deformation associated with intrusions. In a more regional context, the BRVF is situated in a complex tectonic setting that may influence the locations of these intrusions. The regional gravity anomaly and model are explained by thickening of a dense quartzite by thrust faulting (Figure 11c). Such regional density contrasts in the crust are interpreted to influence magma ascent elsewhere (Deng et al., 2017), possibly explained by changes in stress trajectories associated with the differential loads caused by these broad lithologic variations (Connor et al., 2000;Rivalta et al., 2019).

Implications for Volcanic Hazards and Geothermal Exploration
The two anomalies may indicate silicic intrusions occurred at two different times, as indicated by the differing ages of BRVF silicic domes. The CDF alignment erupted approximately 58 ka and the Sheep Island dome, forming an island on the west side of the reservoir, erupted approximately 1.5 Ma (McCurry & Welhan, 2012). This difference in dome ages is consistent with at least two episodes of intrusion. Observations of recent high-level silicic intrusions and eruptions indicate that activity frequently involves a complex series of events (Castro et al., 2016;Jay et al., 2014;Miller et al., 2017;Shaffer et al., 2010). If the intrusions in the BRVF formed coeval with the effusion of the domes, similar to the high-level intrusion at Cordón Caulle (Castro et al., 2016), then it is likely that the northern intrusion was emplaced, in a separate event, prior to the southern intrusion.
The multiple vents of varying ages, the two gravity anomalies and the spatial association with the basaltic volcanic field all indicate that the possibility of future intrusions and dome eruptions should be assessed. Potential for future silicic eruptions in dominantly basaltic volcanic fields changes the way volcanic hazards need to be estimated Duffield et al., 1980;Ewert et al., 2005;Jónasson, 2007;Kósik et al., 2020;Riggs et al., 2019). The CDF events preserve evidence of explosive volcanism, but are comparable or smaller in volume than nearby and more abundant basaltic eruptions. The interpretation of two gravity anomalies as being caused by large-volume and shallow silicic intrusions changes the hazard, since it indicates these eruptive episodes could have evolved into much larger magnitude and intense eruptions with widespread effects. Even as intrusions, deformation appears to be associated with the emplacement of these shallow bodies, and is of much larger amplitude than identified in most basaltic volcanic fields.
Such intrusions and their associated silicic eruptive vents are widespread. Other examples include large-volume exogeneous and endogeous silicic domes erupted on the Eastern Snake River Plain, the Buckskin Dome and Ferry Butte south of the town of Blackfoot and Yandell Mountain southeast of Blackfoot (Figure 1). The CDF domes and tuff rings are small-volume compared to these features (0.46 km 3 ), but the approximately 100 km 3 of the BRVF intrusions is large compared to these other features. For this volume, an intrusive to extrusive ratio for silicic volcanism is 217:1, but recognizing the range of reasonable volumes from the tradeoff curve (Figure 9) gives an intrusive to extrusive ratio can be between 109:1 and 261:1. While the modeled intrusions are high-volume compared with the mapped eruptive products, we note they are less than one-tenth the volume of the largest caldera eruptions and their intrusive magmas (Gregg et al., 2012;Takarada & Hoshizumi, 2020). Eruption magnitudes are classified using orders of magnitude change in volume. In this context, although uncertainty in the volumes of the intrusions is high because of uncertainty in the density contrast, the volume range is consistent with moderately large volume explosive eruptions.
The ages of the eruptions within the BRVF (∼60 ka and 1− 1.5 Ma) suggest a return period of approximately 1 million years. This is a low hazard rate but it also has a high uncertainty with only two constraining events. For comparison, the domes of the ESRP span in from 1.4 − 0.309 Ma and yield a return period of 270 ka (Kuntz et al., 2003). It is possible that the BRVF could have a similar return period for eruptions as the ESRP, considering that the volcanism is chemically congruent (McCurry & Welhan, 2012), but has yet to experience enough volcanism to reflect that similarity. The current geochemical model for the BRVF includes a deeper magma storage system at ∼ 14 km depth (McCurry et al., 2015). Our model is consistent with the conceptual model of a deep magma source. Figure 11c depicts our model of the upper 9 km of the crust spanning from Gem Valley to the northeastern extent of the BRVF. As discussed in previous sections, the gravity high on the west side of the BRVF correlates with the shallowing of Precambrian/Proterozoic quartzites, and the two gravity anomalies in the CDF separated by a saddle are related to shallow silicic intrusions. The wavelengths of the gravity anomalies are short compared to anomalies that would be produced by magma at 14 km depth, and therefore are not related to a deep or midcrustal source. Welhan et al. (2014) investigated heat flow anomalies from the surrounding region to help assess potential geothermal resources and show that the heat flow directly above the modeled intrusions is low while the heat flow to the northwest near the trace of the Meade Thrust is much higher. Given that modeled intrusions are shallow and thin, and were likely emplaced during or before the time of emplacement of China Hat (∼58 ka), they may have completely cooled.

1.
A new gravity survey of the BRVF reveals two negative gravity anomalies underlying and adjacent to late Pleistocene silicic domes and tuff rings. These anomalies, after detrending, have amplitudes up to −16 mgal and ellipsoidal shape, elongated NW. 2. The anomalies are modeled as two shallow silicic intrusions. In map dimensions, each is approximately 9 × 4.5 km. Given the uncertainty in density of the intrusions, their combined volume is estimated to be in the range of 50 − 120 km 3 . Calculated using density contrast of −400 kg m −3 , the northern intrusion has volume of approximately 60 km 3 and the southern intrusion has volume of approximately 40 km 3 . 3. Significant deformation appears to have accompanied the emplacement of these intrusions. NNW-trending fault sets bound the intrusions, with the largest displacement (50 m) observed on any faults in the BRVF immediately adjacent to the southern intrusion. The gravity anomalies are overlain by ENE-trending faults, which may have formed during emplacement and possibly deflation. It is possible that the ascending magma exploited faults in the BRVF and their ascent was influenced by crustal scale structures associated with thrust faults. 4. At least one and likely two episodes of large-volume and shallow intrusion has occurred in the bimodal BRVF. Had these magmas not stalled in the shallowest crust, they would have produced moderately large magnitude eruptions that would have affected broad areas. We suggest identification and quantification of shallow intrusions may help better quantify volcanic hazards in bimodal volcanic fields. Given the tradeoff between density contrast and volume, the intrusive to extrusive volume ratio for silicic volcanism in the CDF is between 109:1 and 261:1.

Data Availability Statement
Datasets for this research are available at Hastings et al. (2021).