The Harney Basin is a relatively flat-lying depression in the northeast corner of the enigmatic High Lava Plains volcanic province in eastern Oregon. A thick blanket of volcanics including flood basalts, rhyolites, tuffs, ash flows, and distinct eruptive centers covers the basin, making it very difficult to study the upper-crustal features. In addition to a portion of the High Lava Plains active source seismic data in the Harney Basin area, we employed geologic, gravity, magnetic, digital elevation, and other geospatial data in our integrated study. We generated an upper-crustal 3D seismic tomographic model of the Harney Basin using a sparse grid of 2D seismic lines and constructed an integrated geophysical model of the upper-crustal structure, which reveals that the basin reaches as deep as 6 km in its central area. The tomographic inversion also detected some unusually high-velocity (>6.5 km/s) bodies in the upper crust near the central basin area. The presence of several ash-flow tuffs and voluminous rhyolites in the Harney Basin region indicates that the sources of these materials are nearby. We observe two major caldera-shaped features within the basin, which we interpret to be likely candidates for the source of some of these tuffs. These potential calderas are associated with low seismic velocities, low gravity anomaly values, and topographic depressions. We interpret the extent and evolution of these potential calderas based on our integrated analysis.


The High Lava Plains (HLP) of Oregon (Fig. 1) are a major component of a large igneous province in the Pacific Northwest that was formed due to complex intra-plate volcanism during the Cenozoic (e.g., Carlson and Hart, 1987; Camp et al., 2003; Streck and Grunder, 2008). The Steens Basalts that cover much of the HLP and the Columbia River Basalts erupted at ca. 16–17 Ma, which was followed by the extensive bimodal rhyolitic and basaltic volcanism in the area (Hooper, 1997). The genesis of the younger HLP time-transgressive track of volcanism (Fig. 1) and its possible relationships with the Snake River Plain (SRP) volcanic track are highly debated. Among many of the hypotheses, backarc extension is postulated as a dominant force (e.g., Eaton, 1984; Carlson and Hart, 1987; Ford et al., 2013).

The HLP area was affected by a variety of tectonic processes including Farallon plate subduction, backarc volcanism, ignimbrite flare ups (Lipman et al., 1971), strike-slip deformation, Basin and Range extension, and terrane accretion in a comparatively short period of time. In the HLP, the time-transgressive Cenozoic volcanic activity covered the area with lava flows and pyroclastic material (e.g., Camp and Ross, 2004; Brueseke et al., 2007; Streck and Grunder, 2008), making it difficult to identify older crustal structures and decipher tectonic events (Anderson, 1989). Thus, the HLP was the focus of a large project funded by the National Science Foundation, which involved a wide range of geophysical and geological studies that helped resolve the crustal structure of the crust and upper mantle (e.g., Roth et al., 2008; West et al., 2009; Eagar et al., 2010, 2011; Druken et al., 2011; Cox et al., 2013).

The Harney Basin, a major physiographic feature in eastern HLP (Fig. 1), has been interpreted as a large structural depression possibly created due to caldera collapse (e.g., Greene, 1973; Walker, 1979; Walker and Nolf, 1981; Streck and Grunder, 1995, 2008). In the Late Cenozoic, the depression was filled with ash-flow tuffs, basalt flows, rhyolites, andesites, basaltic breccias, detrital alluvium, fanglomerate, playa deposits, and lacustrine sediments derived from the widespread igneous rocks (e.g., Piper et al., 1939; Baldwin, 1976; Walker, 1979; Russell, 1984). One element of the HLP Project was to target the Harney Basin with 3D seismic coverage as part of the active source experiment. Gravity data were also collected in the area of the 3D coverage. The purpose of this paper is to present the results of an integrated analysis of the seismic and gravity data that targeted the upper crust of Harney basin. Our approximation of the physiographic boundary of the Harney Basin (Fig. 2) was drawn based on a 10 m resolution digital elevation map of the area and is used in the figures that follow.


The accreted Blue Mountains terrain and the Columbia River Basalt Plateau lie north of the HLP; whereas, to the south the HLP is bounded by the northern Basin and Range Province. The eastern extent of the HLP is truncated at the Owyhee Uplands west of the 0.706 Sr isotope line that separates the North American craton to the east from the Mesozoic accreted terrains (Armstrong et al., 1975), and the Cascade Range lies to the west of the HLP. Other major geological provinces in the area include the Oregon-Idaho grabens, Brothers fault zone, and the Newberry volcanic field (Fig. 1).

In a regional context, the HLP track of bimodal, time-progressive volcanism started near the McDermitt volcanic field (Fig. 1) in southeast Oregon at ca. 16 Ma (Pierce and Morgan, 1992), and it is probably a mirror image of northeast-trending Yellowstone–Snake River Plain volcanic centers. It continued to migrate in a WNW direction with its youngest volcanic center at the recently active Newberry Volcano (e.g., Walker, 1974; Christiansen and McKee, 1978; Smith and Luedke, 1984; Jordan et al., 2004; Meigs et al., 2009). The Brothers fault zone (Fig. 1) is a major structural feature approximately parallel with the HLP track and extends through the Harney Basin. It has been interpreted as one of the major conduits for bands of rhyolitic volcanic deposits in the HLP (Christiansen and Yeats, 1992).


The Harney Basin is one of the main features of the HLP (Fig. 1). It is bounded to the north and south by the Blue Mountains and Steens Mountain, respectively. The eastern side is bounded by the Owyhee Plateau, but the western boundary is ill defined. The basin is mostly covered with Quaternary sediments and seasonal lakes with basalt flows, tuffs, volcanic ash, and volcanic centers.

Three voluminous tuff deposits have been mapped in the Harney Basin region. The oldest (7.1 Ma), Rattlesnake Tuff (RST), is the largest and covers an area of ∼9000 km2 (Streck and Grunder, 2008). When reconstructed, the areal coverage is estimated to be ∼40,000 km2 with a volume of 280 km3 of Dense Rock Equivalent (DRE). Based on the degree of welding, pumice size, and thickness, its source has been inferred to be in the western Harney Basin (Streck and Grunder, 1995, 2008). The second oldest (8.4 Ma) is the Prater Creek Tuff (PCT) and had a DRE volume of ∼100–150 km3 (Greene, 1973; Ford et al., 2013). The youngest (9.7 Ma) is the Devine Creek Tuff (DCT), which covered an area of ∼19,000 km2. Its DRE volume is estimated to be ∼195–250 km3 (Greene, 1973), and its source is interpreted to be beneath the eastern Harney Basin based on the thickness and distribution of the volcanics (Fig. 2). There are numerous small-scale ash-flow deposits in the area with ages ranging from 12 Ma to Holocene.


The HLP active source seismic experiment in 2008 (Fig. 3) was conducted to image the crust and upper mantle structures beneath the HLP and Steens Mountain regions, and it was designed to tie to the previous active source experiments in the region (Catchings and Mooney, 1988 to the west and Lerch et al., 2007, 2008 to the south). The survey consisted of two major ∼400-km-long lines across the HLP and some short linking lines around the Harney Basin that provided a modest 3D seismic survey of the basin. The experiment employed 15 one-ton explosive sources and 2612 Texan seismic recorders with 4.5 Hz geophones. The average spacing between instruments was ∼800 m with an average source spacing of ∼50 km.

In this study, we used six shots around the Harney Basin with five lines of receivers (Figs. 3A and 3B and Supplemental Fig. 11). Raw gathers (super-shot gathers) that contain all recordings for individual shots were generated and consisted of seismograms from 1148 receivers (Supplemental Fig. 22). In order to create a 3D seismic-velocity model of the upper-crustal features beneath the Harney Basin, we subsampled the raw gathers and generated 26 new gathers in the area of interest (Figs. 3A and 3B), and the details are listed in Table 1. We analyzed these seismic data using the fast marching tomographic method (FMTOMO) for 3D tomographic inversions, which has been successfully employed in crustal-scale studies (e.g., Rawlinson and Urvoy, 2006; Rawlinson, 2007; Brikke, 2010; Rawlinson et al., 2010; Rockett, 2011). The major seismic result was a 3D tomographic model of the upper-crustal features.

In addition, gravity, magnetic, geologic, and geospatial data were also employed in our integrated analysis. We generated multiple gravity and magnetic maps as well as integrated gravity models across the Harney Basin. The results from the seismic tomography, gravity models, geospatial data, surface geology, and topography were used to interpret the basin structure and evolution.

The 2008 HLP active seismic data set has been analyzed previously (e.g., Okure, 2009; Cox et al., 2013) and interpreted in the form of 2D seismic lines. Okure (2009) used the HLP seismic data around the Harney Basin area and tentatively located the basin boundary based on the travel-time delays observed in the seismic data. Cox et al. (2013) also used the HLP active seismic data for a crustal-scale study and interpreted the two long 2D lines across the HLP region. A 5–7-km-thick upper layer of sediments and volcanics was detected in the HLP region. In this study, underplating in the lower crust and magmatic modification within the crust were observed beneath the Harney Basin (Cox et al., 2013).

Seismic Data Processing and Preparation for Tomographic Inversion

To enhance the signal-to-noise ratio for the seismic data, spiking deconvolution, trace DC removal, trace equalization, and a 1-4–12-15 Hz Ormsby bandpass filter was applied. Afterwards, we applied crooked line geometry, sorted the receivers and shots, and generated 26 different shot gathers (Figs. 3A and 3B and Table 1). The processing flow is shown in Supplemental Figure 33. We chose the source-receiver pairs to establish the best possible 3D coverage of the upper crust of the basin with the data available. The idealized ray path diagram for the chosen source-receiver pairs is shown in Figure 3B.

An open source code (ZP, http://www.soest.hawaii.edu/users/bzelt/zp/zp.html) was employed for reading, plotting, and picking these SEGY gathers. The first arrival zero crossings before peaks were picked and exported with the picking uncertainty values that are an output of the ZP software (Zelt and Smith, 1992). The details of these data are summarized in Table 1. Examples of the SEGY gathers before filtering are shown in Figures 4A and 4C, whereas processed and picked corresponding gathers are shown in Figures 4B and 4D.

Next, we used the fast marching tomographic modeling package (FMTOMO), to perform the 3D travel-time tomography and to generate seismic-velocity models in and around the Harney Basin area. The tomographic inversion is based on an iterative nonlinear subspace inversion scheme, which is capable of adjusting model parameters in order to satisfy the observed data (De Kool et al., 2006). FMTOMO requires information about source and receiver latitudes, longitudes, and elevations, as well as travel-time picks and associated uncertainty values for the inversion. Initially, we used a simple two-interface velocity model containing 14 nodes with continuous velocity variations between these nodes (Table 2). The model space was set up between 42.5° and 44°N latitude, 118° and 120.25°W longitude, and 1.9 km above mean sea level (bsl) to 15 km depth below sea level (bsl).

A checkerboard test for the model was performed to check the solution robustness of the tomographic inversion. Figure 5 shows the results from a resolution test where we compared the initial and final velocity models before and after inversion. This test shows that the solution is robust and reliable in and around most of the Harney Basin area up to a depth of ∼8 km bsl.

Seismic Tomography Inversion

After the successful checkerboard test, we ran the forward and inverse model in FMTOMO. Stable calculations require an optimized ratio between the grid cell size and ray-path density. In this study, a propagation grid with cell size of 2.6 × 2.3 × 0.55 km and a velocity grid with cell size of 5.2 × 4.6 × 0.85 km were found to be optimal. The details of the tomographic inversion process applied are discussed by Khatiwada (2013) and are also schematically shown in a flowchart (Supplemental Fig. 3 [see footnote 3]). The process of forward and inverse modeling was repeated until satisfactory convergence values were obtained. Out of seven iterative inversions that we ran, the sixth iteration achieved the best convergence (Table 3). During the inversion, root mean square (RMS) travel-time residuals decreased by ∼72%, whereas the variance of the misfit and χ2-misfit values decreased by ∼84% and ∼90%, respectively (Table 3).

We used the Generic Mapping Tools (GMT) software package (http://gmt.soest.hawaii.edu/) to extract various 3D velocity slices and Python programs to visualize these results. Some example velocity slices generated are shown in Figure 6. Figures 6A and 6B are the initial and inverted velocity slices at 3 km depth bsl. Figures 6C and 6D are the initial and inverted velocity models for a longitudinal slice along 119°W; these figures show a significant amount of upper-crustal velocity variation in the Harney Basin area. Similarly, in Figures 6E and 6F, we show the initial and inverted velocity models along 43.60°N latitude.

Seismic Tomography Results

The seismic tomography results successfully imaged the upper-crustal features of the area to ∼8 km depth bsl. Figure 7 shows a series of depth slices ranging from 0 km to 7 km bsl. The velocity normally ranges from ∼2 km/s to ∼7 km/s. Figure 8 shows a series of longitudinal slices from the inverted velocity model. The upper layer basin fill has a velocity of ∼3–3.5 km/s and is ∼2 km thick in the central area, whereas the bottom of the basin is ∼6 km deep in this area. The basin shallows toward the west. On the east, the longitudinal slice at 118.60°W shows that the basin deepens up to ∼5 km bsl.

A series of latitudinal slices are shown in Figure 9. The basin-fill thickness reaches to ∼5 km bsl in the central area and is shallower to the south. The first two layers (yellow and green) represent sedimentary fill and widespread tuffs deposits. The third layer (blue) from the top is most probably the volcanic units (basalt, rhyolites, andesites, and tuffs) that form the base of the basin. Near the central part of Harney Basin (119°W longitude and 43.3°N latitude), we observed an anomalously high velocity (>6.5 km/s) body at ∼5.5 km depth bsl. It lies to the northeast of Diamond Craters area and is about 25 × 20 × 2 km3 in size. Based on the lab experiments, Song et al. (1997) suggested a similar velocity for high-alumina olivine tholeiite (HAOT) basalt at ∼5 km bsl. Thus, we interpret this feature as an intrusive body associated with the low-K, HAOT basalts (e.g., Hart et al., 1984; Hart, 1985).


In addition to the active source seismic data, we used gravity and magnetic data as part of our integrated geophysical interpretation. The gravity and magnetic data for the whole continental United States are freely available in the Pan American Center for Earth and Environmental Studies (PACES) Web sites (http://research.utep.edu/) (Keller et al., 2006). State-based gravity and magnetic data can be found in the U.S. Geological Survey’s database (http://pubs.usgs.gov/ds/355/). We also used digital elevation models (DEMs), geological maps, fault databases, aerial photos, and tectonic maps of the area.

Gravity Data Collection and Preparation

In addition to the existing databases, ∼1000 new gravity data were acquired in and around the Harney Basin area in 2012 with an average station spacing of ∼1.5 km. Similarly, ∼270 gravity data were collected in 2008 in the same area. We used a differential global positioning system (GPS) to obtain the spatial coordinates and elevations of the gravity stations with horizontal and vertical submeter accuracy of ∼30 cm and ∼10 cm, respectively.

The GPS data were processed with the help of Online Positioning Users Service (OPUS), which is an online GPS data processing Web site affiliated with the National Geodetic Survey (https://www.ngs.noaa.gov/OPUS/). All of the data collected were tied to the international gravity base station at Burns, Oregon, to obtain the absolute gravity values.

The gravity reduction spreadsheet of Holom and Oldow (2007) was used to obtain the standard free air and Bouguer anomaly values (2.67 g/cm3 density). The complete Bouguer anomaly (CBA) values were obtained after applying terrain correction to the data. Afterwards, these data were merged with ∼5400 gravity values from the PACES database in the area.

Gravity Anomaly Mapping

We used a 2 km grid spacing to generate a CBA gravity map of the area (Fig. 10A). The anomaly values across the region vary by ∼90 mGal. Major tectonic units and their gravity signatures in the area are easily identified (Fig. 10A). Later, we confined our study area around the physiographic Harney Basin (Figs. 10B, 10C, and 10D).

A series of wavelength and edge-detecting filters were applied in the frequency domain using fast Fourier transform (FFT) techniques. Figure 10B shows the CBA map of the Harney Basin area. The CBA anomaly varies by about ∼50 mGal with lowest values in the Steens Mountain area. An oval-shaped gravity low is observed in the northeast corner of the basin. The Diamond Craters area is also located in a gravity low. The tilt derivative is an edge detector algorithm that was applied to help to enhance the edges of some subtle features as well as the edges of major features such as basins, possible craters, and calderas (Fig. 10C). For example, a circular-shaped feature, possibly a buried crater or caldera rim, with ∼8 km diameter is clearly visible just west of Steens Mountain (Fig. 10C). This previously unnoticed feature could be one of the eruptive centers for the Steens Basalt. We also applied an upward continuation filter to the CBA grid to remove regional anomalies and generated a 15 km residual CBA map by simply subtracting the upward continued grid from the CBA grid (Fig. 10D). The black lines on this figure are the mapped Quaternary surface faults.

Magnetic Anomaly Mapping

As a part of our integrated geophysical interpretation, we also used the U.S. Geological Survey (USGS) state magnetic databases for Oregon and Idaho. The details of the data acquisition, survey parameters, and processing applied to these grids are available in the USGS data repository (http://pubs.usgs.gov/ds/355/). We merged these two grids and produced a reduced-to-pole total magnetic intensity (TMI) map of the area. However, due to the extensive and diverse volcanic activities in the region, the resulting map was too complex to interpret with confidence (Supplemental Fig. 44). The shapes of some anomalies indicate the presence of basins, but we did not significantly depend on the magnetic data in our interpretation.


We modeled three gridded gravity profiles across the Harney Basin, namely X–X′, Y–Y′, and Z–Z′ as shown by dashed white lines in Figure 10D. The seismic-velocity model generated by tomographic inversion along these profiles, surface geology, and the digital elevation model (DEM) data were used as the key constraints for these gravity models. To find the density (ρ) of the modeled units, we used the relation between P-wave velocity (Vp) and density postulated by the Nafe-Drake curve and given by Brocher’s (2005) relation: 

All of the profiles were modeled to a depth of 15 km. Figure 11A is a gravity model along X–X′ (43.25°N latitude). The two uppermost layers with densities ranging from 2.43 to 2.54 g/cm3 are the Quaternary sediments with occasional ash flow, tuffs, and flood basalt. These layers reach as deep as ∼3 km bsl. The third and the fourth layers with densities of 2.64 and 2.69 g/cm3, respectively, are mostly mixed lithologies of flood basalt, rhyolites, and sedimentary units. The fourth layer is modeled as deep as 8 km bsl and is interpreted to be the crystalline basement.

Figure 11B is a gravity model along Y–Y′ (43.60°N latitude). The eastern portion of the model is a basinal area covered by a thick upper layer (yellow) with a density of 2.43 g/cm3 as well as a thinner low-density layer (green). The total thickness of these layers reaches a maximum of ∼4 km bsl, and as discussed later, we interpret this area as a possible caldera. Beneath the Harney Lake area, at ∼7 km depth bsl, we included a body with a density of 2.79 g/cm3 that corresponds with a high-velocity body in approximately the same location. As discussed above, we interpret it as a residual mafic body that possibly consists of HAOT-type basalt. A horst block related to the Brothers fault zone is also included in this model based on geologic data.

The gravity model along line Z–Z′ (118.8°W longitude, Fig. 11C) crosses two basinal areas (Fig. 10D), which are as deep as ∼6 km bsl. We interpret these features to be possible calderas that are separated by a gravity high near Harney and Malheur Lakes (G1, Fig. 10D). This gravity high is caused by a high-density (2.78 g/cm3) mafic body at the depth of ∼4 km bsl.

These gravity models show the complexity of the subsurface of the Harney Basin and match approximately with the seismic models. Also, the magnetic profile in Figure 11C generally matches with the modeled subsurface lithologies specifically above the interpreted calderas.


Based on these geophysical results and geological observations, we present additional evidence in support of the caldera interpretations. The presence of the late Miocene (10–5 Ma) Rattlesnake, Devine Canyon, and Prater Creek Tuffs, their reconstructed volumes as well as other numerous smaller tuff deposits in the area indicate that their sources may be in the Harney Basin area. The circular to oval topographic lows and the abundance of tuffs, ash flows, scoria, and rhyolites make Harney Basin a suitable place to contain the sources of these deposits. Previous researchers have also interpreted that the sources of these tuffs lie within the Harney Basin (e.g., Parker, 1974; Walker, 1974, 1979; Macleod et al., 1975; Walker and Nolf, 1981; Streck and Grunder, 1995, 2008, 2012; Ford et al., 2013). The source locations of the tuffs were interpreted based on the degree of welding, the changing thickness of the tuffs, changing pumice size, and the shape of the tuff fragments. After a large volume of material erupted, the bimodal magmatism and sedimentation continued in the area for a long time and covered most of the previous igneous features. By combining this geologic information and the results from our gravity and seismic tomographic analyses, we suggest the presence of calderas in the area and discuss this possibility further below.

Integrated Geophysical Interpretation

The gravity anomalies represent density variations of subsurface bodies and, thus, should have a direct relation with the P-wave velocities. In order to understand the relationships between the observed gravity anomalies and the velocity model, we generated 3D block diagrams. Figure 12A illustrates the residual Bouguer anomaly map on the top and a seismic velocity slice through 119.0°W longitude in the cross-sectional view. The interpreted physiographic boundary of the Harney Basin is shown by the yellow line, and its southern boundary is associated with an upwarp of high-velocity values that indicates a structural boundary. We observed a high-velocity (>6.5 km/s) body near the center of the vertical slice at about a depth of ∼5 km bsl. This feature coincides with a high-gravity anomaly (G1 in Fig. 12A) in the Harney Lake area. The northern part of the basin contains an oval area of low-gravity values (Fig. 12B), which coincides with the thicker postcaldera basin-fill deposits shown by the dashed black lines in Figure 12A. The dashed white line is the possible base of the basin that formed the base of the caldera.

Figure 12B is another block diagram with the residual Bouguer anomaly map on the top and a seismic velocity slice through 43.10°N latitude in the cross-sectional view. We observed two elliptical-shaped gravity lows whose shape resembles that of calderas in other volcanic provinces. Hereafter, we refer to these calderas as the northern caldera and the southern caldera (Fig. 12B). We also observed a high-velocity anomaly at ∼4 km depth bsl beneath the Diamond Craters area indicating the presence of a mafic intrusion.

For a detailed integrated geophysical interpretation within the Harney Basin, we compared the results from the gravity and seismic data, a DEM, and magnetic data as shown in Figure 13. Figures 13A–13D are a residual Bouguer anomaly map, a 15 km upward continued residual RTP TMI map, a 10 m resolution DEM map, and a seismic-velocity slice at 4 km depth bsl, respectively. The black dots in Figure 13A are the gravity stations, and the white line in Figures 13A–13C is the interpreted physiographic boundary of the Harney Basin (Fig. 2). The residual gravity anomaly map in Figure 13A contains two elliptical gravity lows bounded by highs that are interpreted as the northern and the southern potential calderas as shown by dashed yellow and white lines, respectively. These areas also possess anomalously low seismic velocity at shallow levels and low magnetic anomalies (Figs. 13B, 13C, and 13D). We interpret the NW-SE–trending linear gravity high in (A), a magnetic high in (B), and a linear velocity high in (D) shown by a pair of black arrows as mafic dikes. The trend of these features is close to that of the Brothers fault zone. We interpret this feature as one of the possible sources for magmatic activity in the Harney Basin area. In order to further investigate these possible caldera features, we divided the basin into two areas as discussed below.

Diamond Craters and Southern Caldera

The Diamond Craters is a complex of basaltic lava flows of late Pleistocene to Holocene age that lies in the southern part of the Harney Basin and covers ∼60 km2 area (e.g., Piper et al., 1939; Peterson and Groh, 1964). We used gravity, magnetic, seismic, and DEM data in our detailed integrated geophysical interpretation of the area and show the results in Figure 14. Figure 14A is a residual Bouguer anomaly map of the area draped on a 10 m resolution DEM map. The Diamond Craters location is shown by a red dashed line. To the west of the Diamond Craters, we observe an elliptical-shaped gravity low shown by the dashed yellow line. This interpreted caldera has diameters of ∼10 km by ∼5 km. We plotted three seismic-velocity slices (Figs. 14B–14D): a horizontal slice at 4 km bsl (B), a longitudinal vertical slice through 118.90°W (C), and a latitudinal vertical slice through 43.10°N (D). In each of these images, comparatively low velocity material extends as deep as ∼4.5 km bsl suggesting the presence of a caldera. This interpretation is also supported by the topographic features, surface geology, and magnetic data (Figs. 10 and 12B). We propose this possible caldera to be the source of 9.7 Ma old Devine Canyon Tuff. This location is close to one proposed by Greene (1973) and Walker (1974, 1979). Beneath the proposed caldera, we observed an ∼1.5-km-thick, high-velocity body at a depth of ∼6 km bsl.

Based on these results, we propose a three-stage evolution of the Diamond Craters and southern caldera area. At ca. 10–8 Ma, the area was covered with active felsic volcanism (e.g., Walker, 1979; Christiansen and Yeats, 1992; Jordan, 2002). Co-eruptive subsidence occurred, and the southern caldera formed. The caldera was later filled by basin sediments, tuffs, and ash flows. Later, due to the faulting in the Brothers fault zone, the olivine-rich basaltic magma (Russell and Nicholls, 1987) rose along these fissures. During the Pleistocene, this magma made its way to the surface in the eastern portion of the caldera, and the Diamond Craters were formed. Some of the residual magma might have solidified in the upper crust at about ∼6 km depth bsl in the nearby area. Meanwhile the sedimentation and volcanism continued to form the present-day landscape in the area, probably leaving a buried caldera beneath the Diamond Craters area.

The Northern Caldera

The northern proposed caldera (∼20 km by ∼30 km) lies in the northeast corner of the Harney Basin and includes the northern half of Malheur Lake (Fig. 15). It is bounded to the north and east by topographic highs. The boundary of the northeast quadrant of the proposed caldera is particularly well defined by a topographic discontinuity. The proposed caldera boundary is delineated with the dashed yellow line. The caldera is elliptical in shape and contains two distinct depressions shown by C1 and C2 (Fig. 15A). The west and south sides are relatively flat. The area is covered by fluvial, fluvio-lacustrine, and volcanic sediments.

We analyzed the gravity, magnetic, DEM, and seismic data in our detailed study. Figure 15A shows the residual Bouguer anomaly draped over a 10 m DEM map. Figures 15B, 15C, and 15D are seismic-velocity slices at 4 km depth from bsl, a longitudinal vertical slice at 118.90°W longitude, and a latitudinal vertical slice at 43.55°N latitude. The low-velocity material has a thickness of ∼4–5 km (Figs. 15C and 15D), and we interpret this as the sediments and volcanic material that filled the caldera after it was formed. Thus, the surface geology, topographic features, and geophysical data suggest the presence of caldera in the area, and in addition, Greene (1973), Walker (1974), and Ford et al. (2013) suggested the source of the 8.4 Ma old Prater Creek Tuffs was present in this same area.

One can argue for alternative interpretations to the presence of caldera such as: thick basin fill, the presence of a silicic batholith, or both, which could produce the low-gravity values and seismic velocities as well as magnetic signatures. Although sparse (Fig. 2), the drilling data suggest the sediment fill in the basin is not more than a few hundreds of meters thick. We do not see any evidence of Basin and Range structures such as fault-bounded grabens and significant extension. There is no evidence for present or past large paleofluvial systems that could fill the Harney Basin with fluvial sediments. Finally, the presence of unusually high-density and seismic-velocity material in the shallow crust (Figs. 7–13), lack of extensive granitic rocks in the vicinity, and the presence of materials that produce higher magnetic anomalies in the caldera area indicate the absence of granitic plutons.

The presence of abundant rhyolite raises the question of its source. Based on the research on the Rattlesnake Tuff, these rhyolites and tuffs are rich in iron and silica content, have high parent magma temperatures (>880 °C), and have high Ba/Sr ratios (e.g., Streck and Grunder, 1995, 2008, 2012; Ford et al., 2013). These characteristics indicate that the sources of these deposits were mafic, and fractional crystallization created the large volume of tuffs and rhyolites, thus ruling out the possibility of the presence of silicic plutons.


Using a series of 2D seismic lines to generate a 3D seismic tomographic model of upper-crustal structure in the Harney Basin area produced meaningful results that were cost effective. The tomography model, gravity and magnetic maps, integrated gravity models, surface geology, and other geospatial data were combined to interpret the major structural components of the upper crust in the Harney Basin area. The newly acquired gravity data revealed some previously unseen features in the study area. Our integrated geophysical approach helped to image the detailed upper-crustal features of the Harney Basin area down to a depth of ∼8 km bsl. Our new integrated geophysical results helped us to identify two potential caldera complexes within the Harney Basin, namely the southern and northern calderas.

Previous researchers also suggested the presence of calderas (Fig. 2) in the area (e.g., Greene, 1973; Walker 1979; Walker and Nolf, 1981; Streck and Grunder, 1995, 2008; Ford et al., 2013) based on geological data. The presence of tuff deposits and numerous ash flows within the Harney Basin is additional evidence to support our interpretation. The complexity of the caldera shape and features suggests that there is no single end-member mechanism for the caldera-forming processes. The sizes of these potential calderas are comparable to the Valles Caldera in New Mexico (Ankeny et al., 1986; Olsen et al., 1986). However due to the extensive crustal modification and magmatic activities in the HLP during late Tertiary and Quaternary, the evolution of these calderas is much more complex in comparison to the Valles Caldera.

To strengthen the proposed geophysical interpretation and further understand the complex evolution of the upper-crustal features and proposed calderas within the Harney Basin, scientific drilling of these features and a detailed 3D seismic study of the area would be ideal.

We are grateful to Dr. Nick Rawlinson for providing us with useful insight about FMTOMO and interpretation of some results. We thank Prof. Anita Grunder and Dr. Benjamin Drenth for their constructive and encouraging suggestions in regard to our interpretations. We also thank those who helped collect the gravity data in the field and the reviewers of this paper for their thoughtful comments. This study was funded by National Science Foundation Continental Dynamics Program grant EAR-0641515.

1Supplemental Figure 1. Experimental layout for the High Lava Plains (HLP) active seismic project conducted in September 2008. Yellow stars represent source locations, whereas the purple and blacks dots are the receiver locations. Fifteen shots (1 ton each), 2612 Texans recording geophones, and 120 RefTEK130 receivers were deployed in this survey. If you are viewing the PDF of this paper or reading it offline, please visit http://dx.doi.org/10.1130/GES01046.S1 or the full-text article on www.gsapubs.org to view Supplemental Figure 1.
2Supplemental Figure 2. Flowchart used for this integrated geophysical interpretation of the Harney Basin region. The seismic, gravity, and magnetic processing flows are independent, while the interpretations are integrated. If you are viewing the PDF of this paper or reading it offline, please visit http://dx.doi.org/10.1130/GES01046.S2 or the full-text article on www.gsapubs.org to view Supplemental Figure 2.
3Supplemental Figure 3. An example of the raw super-shot gather for shot point 24 from the High Lava Plains seismic experiment. The super-shot gather was generated in such a way that it contains all the receivers (1148) that recorded the particular shot. If you are viewing the PDF of this paper or reading it offline, please visit http://dx.doi.org/10.1130/GES01046.S3 or the full-text article on www.gsapubs.org to view Supplemental Figure 3.
4Supplemental Figure 4. Reduced-to-pole total magnetic intensity map of the study area in Figure A and a map of Harney Basin area in Figure B. The black polygon is the physiographic boundary of the Harney Basin. The black lines in Figure A represent mapped faults in the area. The magnetic anomalies related to Blue Mountain terranes, Western Snake River Plain (WRSP), Owyhee Plateau, and High Lava Plains (HLP) are identified. In Figure B, possible calderas within the northern Harney Basin and the Diamond Craters area are shown. The Steens Dikes (SD) are indicated by a pair of dashed black lines. If you are viewing the PDF of this paper or reading it offline, please visit http://dx.doi.org/10.1130/GES01046.S4 or the full-text article on www.gsapubs.org to view Supplemental Figure 4.