Geophysical study of the San Juan Mountains batholith complex, southwestern Colorado

The topographically elevated and deeply dissected San Juan Mountains of southwestern Colorado, USA (Fig. 1), are principally the erosional remnants of the mainly Oligocene San Juan volcanic field (SJVF) (Fig. 2), the laterally and volumetrically largest subarea of the Southern Rocky Mountain volcanic field (Steven, 1975; Lipman, 2007). The San Juan Mountains are bordered by the topographic and structural San Juan and San Luis Basins, to the south and east, respectively, the Colorado Plateau to the west, and the high Rocky Mountains of central Colorado to the north (Fig. 1).

One of the largest and most pronounced Bouguer gravity anomaly lows over North America is coincident with the San Juan Mountains (Fig. 3). This anomaly has been interpreted as the manifestation of a low-density, upper crustal granitic batholith complex that represents the plutonic roots of the SJVF (Plouff and Pakiser, 1972). Whereas this interpretation remains essentially unchallenged, new gravity data processing techniques, digital elevation data, and constraints from seismic refraction studies (Prodehl and Pakiser, 1980) enable reassessment and improvement of the previous model.

We employed a three-dimensional forward modeling approach to analyze the gravity signal of the region using density constraints from previous studies (Table 1) and topographic relief represented by digital elevation data. These data are expressed in a modified Bouguer-like regional gravity map. The upper crustal density structure is significantly better constrained using this modified gravity map, resulting in modification of both lateral and depth extents of the modeled batholith complex and estimation of its total volume. Our analyses are undertaken within the context of recent interpretations of underlying mantle structure (Roy et al., 2004) and petrologic evolution of the SJVF (Lipman, 2007; Farmer et al., 2008).

The modern San Juan Mountains are underlain by geologic units ranging from Precambrian crystalline rocks to middle to late Cenozoic volcanic rocks of the SJVF that dominate the surface geology (Fig. 2). Structurally high Precambrian rocks crop out to the north and southwest of the SJVF (Gunnison uplift and Needle Mountains, respectively; Fig. 2). The oldest exposed Precambrian rocks in the study area are ca. 1.7 Ga metasedimentary and metavolcanic rocks (included in undifferentiated Precambrian rocks; Fig. 2). Subsequent magmatism ca. 1.7 Ga and ca. 1.4 Ga left behind large plutons and possibly batholiths in the region (units Xg and Yg; Fig. 2) (Tweto, 1979; Bickford and Anderson, 1993; Anderson and Cullers, 1999; Karlstrom et al., 2004; Gonzales and Van Schmus, 2007). These plutons may be volumetrically significant in the subsurface, although their lateral and vertical extents are unknown.

A number of tectonic events since the Precambrian have affected the SJVF region. Cambrian rifting in Oklahoma, reflected as the southern Oklahoma aulacogen, may have extended into and through southern Colorado, including the SJVF region (Larson et al., 1985; Casillas, 2004). Rocks associated with rifting include mafic and ultramafic intrusions that crop out immediately north of the SJVF (Fig. 2). Phanerozoic sedimentary rocks of the Colorado Plateau are along the western margin of the San Juan Mountains, and the region was affected by the Ancestral Rocky Mountain orogeny (Kluth and Coney, 1981; Kluth, 1986), during which the Uncompahgre uplift formed. This uplift extends to the northwest of the SJVF (Casillas, 2004). The Laramide orogeny was responsible for uplift of the Needle dome and Gunnison uplift (Kelley, 1955; Tweto, 1975; Cather, 2004) and formation of the San Juan sag, a northward extension of the San Juan Basin into the eastern SJVF region that is filled with as much as ∼2 km of Laramide-age (ca. 70 to ca. 40 Ma) sedimentary rocks (Gries, 1985; Brister and Chapin, 1994).

Most of the study area is underlain by the SJVF, a 25,000 km² erosional remnant of an Oligocene volcanic province that was active across the southern Rocky Mountains of southern Colorado and northern New Mexico (Steven, 1975; Lipman, 2007). Magmatism that produced the SJVF is interpreted as having resulted from a long-term, complex interaction of mantle source melts and continental crust related to subduction beneath the western U.S. (Lipman et al., 1978; Farmer et al., 2008). The SJVF is predominantly andesitic, but ranges compositionally from basalt to rhyolite (Lipman et al., 1970; Lipman, 2007, and references therein). Early andesitic volcanism was followed by voluminous ignimbrite eruptions of dacite and rhyolite, which resulted in the formation of multiple nested caldera complexes (see following). Mixing of basaltic magmas with crustal rocks is inferred to be the source of the silicic products emplaced into and erupted from the upper crust (Lipman et al., 1978). The minimum volume of basaltic magma required to be emplaced into the crust to produce the SJVF rocks is estimated to be 384,000 km³ (Farmer et al., 2008). There may be a significant amount of mafic residuum left behind in the upper mantle and lower crust, based on analogy with other regions of the southwestern U.S. (Keller et al., 2005) and bounds on possible magma volumes from petrologic modeling (Farmer et al., 2008).

SJVF volcanism began ca. 35 Ma, with eruption of andesitic and dacitic lavas of the Conejos Formation from stratovolcanoes until ca. 30 Ma, when eruption of silicic ash-flow tuffs from numerous calderas began (Lipman et al., 1970; Lipman, 2007, and references therein). The Conejos Formation accounts for about two-thirds of the pre-erosional volume of the volcanic field (unit Tpl; Fig. 2) (Lipman et al., 1970). The younger caldera complexes formed within and around clusters of these stratovolcanoes (Steven and Lipman, 1976). Silicic ignimbrites (unit Taf; Fig. 2) erupted from 18 calderas and range in estimated volume from 25 to 5000 km³ (Lipman, 2000, 2007; Lipman and McIntosh, 2008), and individual eruptions are thought to have been among the most violently explosive events in Earth's history (Mason et al., 2004). The first eruptions, from 33.7 to 29.8 Ma, occurred in the northeastern portion of the SJVF, followed by eruptions from 29.4 to 28.8 Ma at the Platoro and Summitville calderas in the southeast SJVF. Volcanic activity then shifted to the western caldera cluster with formation of five calderas between 28.6 and 27.6 Ma. The massive La Garita caldera erupted at 27.8 Ma, the first of 7 nested calderas to form in the central cluster, ending with the Creede caldera at 26.9 Ma (Lipman, 2000, 2006). Eruptions of andesites similar to the Conejos Formation continued throughout the stage of ignimbrite volcanism, such that these rocks are interbedded with the ignimbrites.

Caldera-centered eruptions of ignimbrites continued until ca. 26 Ma, when volcanism changed to a volumetrically minor assemblage of basalts and rhyolites, related to the formation of the Rio Grande rift (Steven and Lipman, 1976; Lipman, 2000, 2007). The Lake City caldera formed at 23 Ma as part of this latter episode within the preexisting western caldera cluster. The ignimbrites (including those of the Lake City caldera) account for about one-third of the original total volume of the SJVF (Lipman et al., 1970), but have been eroded to the point that they now only represent about one-sixth of the total volume (P. Lipman, 2007, written commun.). The total volume of volcanic rocks erupted in the SJVF, including both the Conejos Formation and ignimbrites, is ∼40,000 km³ (Lipman, 2007, and references therein).

A large-amplitude Bouguer gravity anomaly low occurs over the SJVF (Figs. 3 and 4), and has been interpreted to reflect a low-density, composite granitic batholith in the upper crust (Plouff and Pakiser, 1972). This interpretation is discussed in greater detail in the following, but has been accepted in a broad sense by subsequent studies in the region. Such a batholith complex could be considered the accumulated remnants of magma chambers that sourced caldera-forming eruptions (Lipman, 1984). A more recent interpretation holds that batholith complex growth continued incrementally after caldera subsidence ended, as recorded by ages and compositions of localized, shallow intrusions (Bachmann et al., 2007; Lipman, 2007). A number of granitic intrusions crop out in the SJVF region (unit Tmi; Fig. 2), and are associated with the ignimbrite stage of volcanism. They are as much as millions of years younger than the volcanic units at the corresponding calderas, and tend to be more mafic in composition than the volcanic rocks, especially those that are significantly younger than the ignimbrites (Bachmann et al., 2007; de Silva and Gosnold, 2007; Lipman, 2007). A subvolcanic batholith complex remaining behind today is not necessarily compositionally or volumetrically representative of large magma chambers that existed during or immediately following volcanism. It follows that a batholith complex would not be expected to be compositionally uniform, because it would have been assembled over a long period of time of magmas with differing compositions, with the possibility of relatively mafic lower portions (Bachmann et al., 2007; de Silva and Gosnold, 2007; Lipman, 2007).

Gravity anomalies reflect lateral variations of density, with gravity highs occurring over areas of relatively high rock density, and gravity lows occurring over areas of relatively low density rocks. Regional quality gravity data, with station spacings of 1–10 km, were extracted from the PACES (Pan American Center for Earth and Environmental Studies) database (maintained by the University of Texas at El Paso; http://gis.utep.edu, accessed December 2008) and supplemented with 323 stations recently acquired in the San Juan sag region (R. Gries, 2007, written commun.). The PACES database consists of data collected over decades by many previous workers and was compiled as the result of a major cooperative effort between federal agencies and universities (Keller et al., 2006). For comparison to the nonstandard approach taken in this report, a standard Bouguer anomaly map was computed using conventional techniques (Telford et al., 1990; Blakely, 1995). Corrections included those for predicted gravitational attraction at the elevation and latitude of the gravity stations (free air and theoretical corrections), effects of tabular and homogeneous rocks masses between the stations and sea level (Bouguer slab correction), and effects of topographic masses (terrain corrections). The standard reduction density of 2670 kg/m³, assumed to represent an average upper crustal density (Hinze, 2003), was used for the Bouguer and terrain corrections. Application of these corrections yields complete Bouguer anomalies (Fig. 3).

The wavelengths of gravity anomalies can be related to the depths of the causative sources, because deep and/or broad sources produce longer wavelength anomalies than those caused by shallow and/or narrow sources. One method used to enhance the signature of deep and/or broad sources at the expense of shallow and/or narrow sources is upward continuation, in which the gravity anomalies measured on the ground surface are mathematically transformed to a higher measurement level (Blakely, 1995). Here the effect of the gravity field was calculated at a level of 10 km above the ground, in order to enhance the signature of relatively deep/broad geologic features (Fig. 5).

Seismic refraction data indicate the broad velocity structure of the crust and upper mantle, highlighting variations with depth. Whereas these data sets do not provide detailed images of lateral variations of velocities or of stratigraphic layering, they do give direct estimates of seismic P-wave velocity for geologic situations where velocity increases with depth. A seismic refraction profile was recorded over the eastern SJVF by the U.S. Geological Survey during 1965, with shotpoints on the south and north flanks of the SJVF (Fig. 1) (Prodehl and Pakiser, 1980). By modern standards these data are low quality, with low signal to noise ratio and a small number of imprecisely located receivers, yet they are still useful for studies of gross velocity structure. Seismic velocities typically increase with depth, although in the southern Rocky Mountain region, anomalously low velocities have been observed in the crust and mantle, geographically correlated with granitic batholiths (Schneider and Keller, 1994).

Densities of igneous rocks vary as functions of both composition and texture (Telford et al., 1990). Silicic rocks generally have lower densities than mafic rocks, and for a given composition volcanic rocks typically have lower densities than their plutonic equivalents (Telford et al., 1990). Granitic plutons commonly have lower densities than the country rocks they intrude, producing gravity lows (Bott and Smithson, 1967). A number of Bouguer gravity anomaly lows in the southern Rocky Mountain region have been interpreted to reflect low-density granitic batholiths and plutons, including an Oligocene batholith in central Colorado (Isaacson and Smithson, 1976), Oligocene and Precambrian plutons in the Sangre de Cristo Mountains (Cordell and Keller, 1984; Cordell et al., 1985; Grauch and Keller, 2004; Quezada et al., 2004), and an Oligocene batholith under the Mogollon-Datil volcanic field (Schneider and Keller, 1994).

The Bouguer gravity anomaly low associated with the SJVF has dimensions of ∼100 by 150 km and laterally is contained nearly completely within the boundaries of the SJVF (Figs. 2 and 3). An upper crustal composite Oligocene granitic batholith complex related to the SJVF has been interpreted as the source of the anomaly (Plouff and Pakiser, 1972). A number of lines of evidence support this conclusion (Plouff and Pakiser, 1972; Bachmann et al., 2007). (1) The boundaries of the gravity low correlate well with the mapped extent of calderas (Fig. 2), which would be expected for a batholith complex that is related to surface manifestations of magmatism (Lipman, 1984). (2) Low-density volcanic rocks of the SJVF extend well beyond the area of the gravity low and tend to have small density contrasts with surrounding sediments and sedimentary rocks (Table 1) of the Colorado Plateau, San Juan Basin, and San Luis Basin, meaning that the surficial rocks cannot be the main source of the anomaly. (3) Whereas there is a profound inverse correlation between surface elevation and the complete Bouguer gravity anomaly values (Fig. 3), the gradients at the margins of the gravity low are too sharp to be caused by a source >∼15–20 km depth at most and probably 5 km or less on average (Plouff and Pakiser, 1972). This means that Airy-type isostatic compensation (e.g., thickened crust under mountains ranges) cannot be responsible for a significant portion of the gravity low, unlike anomalies observed over many other mountain ranges (e.g., Blakely, 1995). (4) Precambrian granites exposed in the Needle Mountains and in the Gunnison region may be expected to contribute to the gravity low, as is the case in the southern Sangre de Cristo Mountains (Quezada et al., 2004). However, the Bouguer gravity anomaly low does not extend over the entire Needle Mountains, only over its northern portion (Figs. 2 and 3), and at local scale, Precambrian granites there do not appear to correlate with gravity anomalies. This led Plouff and Pakiser (1972) to interpret that the batholith complex was intruded underneath the northern portion of the Needle Mountains. The gravity low also does not extend over granitic rocks in the Gunnison region. Whereas Precambrian granites may account for some of the gravity low and their presence under the SJVF cannot be ruled out, it seems unlikely from these observations that they are a significant source for the gravity low.

Individual calderas are not associated with gravity anomalies within the larger gravity low, indicating that the densities of rocks that fill the calderas are not significantly different from the surrounding volcanic rocks. The Cochetopa caldera is an exception, with a geographically coincident gravity low (Figs. 2 and 3). The source of the low is unknown, but may be caused by a density contrast with relatively dense Precambrian rocks surrounding much of the caldera margin (Plouff and Pakiser, 1972).

Plouff and Pakiser (1972) placed the margin of the batholith complex at approximately the –300 mGal contour. They used two-dimensional forward modeling (close to profile A–A′ of this study) and an assumed density contrast of –100 kg/m³ to estimate a batholith complex thickness of ∼24 km (Plouff and Pakiser, 1972), but did not estimate its volume.

A regional gravity low, 400–500 km in wavelength and therefore encompassing the gravity low related to the batholith complex, is over the region surrounding the SJVF (Fig. 5). Roy et al. (2004) interpreted the source of this low as a low-density zone in the upper mantle, developed as a product of Tertiary magmatism by extraction of basaltic magmas from the mantle into the lower crust. A solution for the source of the long-wavelength gravity low is an ∼300-km-wide, ∼200-km-thick zone assuming a density contrast of –33 kg/m³ with surrounding mantle rocks, on the basis of a model that corresponds in location to profile B–B′ herein (Roy et al., 2004).

Seismic low-velocity zones are present in the upper crust of the SJVF and Aspen-Gunnison regions. Seismic refraction data collected by the U.S. Geological Survey during the 1960s revealed significant thicknesses of relatively low velocity rocks within the upper 20 km of the crust (Prodehl and Pakiser, 1980). Rapid terminations of upper crustal refractions were interpreted to reflect two stacked and laterally extensive zones (with velocities <6 km/s between rocks with velocities >6 km/s), including a small zone of high-velocity rocks between the low-velocity zones (Prodehl and Pakiser, 1980). Whereas the rocks associated with the low-velocity zones were interpreted to occupy the same upper ∼20 km of the crust occupied by the batholith complex beneath the SJVF (Prodehl and Pakiser, 1980), we are not aware of any literature that has commented on the relationship between the source of the seismic low-velocity zones and the source of the gravity low.

Gravity data processing (or reduction) involves the removal of several predictable effects that influence observed gravity values, with the goal of enhancing the geologic signal to be investigated. As discussed herein, conventional gravity maps (such as that shown in Fig. 3) are derived and interpreted assuming that an upper crustal density of 2670 kg/m³ (Hinze, 2003) is appropriate for the Bouguer slab and terrain corrections over large areas.

Interpretation problems arise from conventional reduction procedures when the densities of rocks above sea level differ from 2670 kg/m³. The Bouguer slab and terrain corrections are of special interest in this study, because the andesites and ignimbrites that compose the high elevations (Figs. 2 and 4) of the SJVF are much less dense than the standard reduction density of 2670 kg/m³, and because the San Juan Mountains contain extremely rugged topography and therefore large-magnitude terrain corrections are required. The Bouguer slab correction is proportional in magnitude to the assumed density and algebraically negative in effect, meaning that calculated Bouguer gravity anomalies will be lower for larger assumed densities. Intermediate and silicic composition volcanic rocks, like those that compose the SJVF, normally have densities significantly <2670 kg/m³. This means that the gravity low over the SJVF (Fig. 3) has a significant contribution from the volcanic rocks, and is not entirely due to the batholith complex. Thus, quantitative interpretation of the anomaly will result in models of the batholith complex that are incorrect, if the entire gravity low is attributed to the complex.

To compute a gravity anomaly map that better reflects the batholith complex, we used a method that combines data reduction (corrections similar to standard Bouguer slab and terrain corrections) and modeling of surface geology into one step, and allowed variable surface densities to be used (Wingrav program, M. Baker and K. Crain, 2003–2005, personal communs.). This method calculates the gravitational effect of mass between two layers at each gravity station, where the layers are defined by digital elevation models. Different calculation methods, following those of Hammer (1974), were used depending on the distance from the gravity station to different masses, where boundaries between different masses are defined in terms of a grid of density values (Fig. 6). The edges of individual masses were treated as vertical boundaries, extending from the top of the upper layer to the top of the bottom layer (Crain, 2006). For masses >500 km away from a gravity station, a row mass approximation (Hammer, 1974) was used, and for distances <500 km to a gravity station, vertical line element approximations (M. Baker, 2003, written commun.) were employed. Groups of gravity stations were combined, or binned, into groups of masses for ease of calculation, with the spatial extent of station bins depending on the distance to masses.

Density information for the model comes from published data and typical values based on lithologies (Popenoe and Steven, 1969; Popenoe and Luedke, 1970; Grauch and Hudson, 1987; Jenkins, 1989; Telford et al., 1990), as well as information from borehole density logs (R. Gries, 2007, written commun.). The most important density assignment was that used for the SJVF, because these volcanic rocks form most of the high topography in the study area. Andesites of the Conejos Formation, as well as andesites that continued to erupt during the ignimbrite phase of volcanism, have typical densities of 2500 kg/m³ (R. Gries, 2007, written commun.) and make up ∼5/6 of the total SJVF volume. Ignimbrites make up 1/6 of the SJVF volume and have typical densities of 2200 kg/m³. A weighted average of 2450 kg/m³ was used to represent the entire SJVF, because the subsurface distribution of the different volcanic rocks is not known well enough to be built into the model. Geologic units of small thickness, such as Rio Grande rift–age volcanic flows, were not included in the model because they are volumetrically insignificant. Density assignments are summarized in Table 1.

For this study, we ideally wanted to calculate the gravitational effect of all masses above the level of the batholith complex, so that this effect could be removed and the resulting gravity anomalies would not include density variations related to surficial rocks. In practice, however, the depth to the top of the batholith complex is not well known, nor is the overlying three-dimensional distribution of densities. Also, adjacent to the SJVF are a number of complex structures, including the San Juan and San Luis Basins and San Juan sag, that we did not attempt to model in three dimensions. Thus, we proceeded by assigning variable densities (Fig. 6) and calculating the gravitational effect of rocks only above an elevation of 2286 m (7500 ft), the approximate lowest surface elevation of the adjacent San Luis Basin. The structures of basins and uplifts of Precambrian basement rocks around the margin of the SJVF were not built into the model. For simplicity we further assumed that all density boundaries modeled above 2286 m extended vertically down to 2286 m, such that no dipping contacts are used. All rocks beneath 2286 m were assigned a density of 2670 kg/m³. The gravitational effect was calculated for this model, the upper surface of which is defined by 90 m Shuttle Radar Topography Mission digital elevation data. This effect was subtracted from the original free air gravity field (corrected for distance above the center of the Earth but not masses of rock), in order to yield a residual map (Fig. 7) that is conceptually similar to the standard complete Bouguer anomaly map (Fig. 3). To summarize, data reduction and three-dimensional modeling have effectively been combined to remove the effect of rocks above 2286 m. Thus, this map gives an improved representation of buried anomaly sources.

The difference (Fig. 8) between the standard complete Bouguer anomaly map from the new map is as much as 10–15 mGal, the greatest difference occurring over the highest elevations of the SJVF region (Figs. 4 and 8). The cumulative error level within the new gravity map is estimated to reach ∼3 mGal based on observations of solutions at single stations thought to be spurious, likely the combined result of using relatively coarse (90 m) digital elevation data and the inherent resolution of the Wingrav program with the chosen calculation parameters. The map is therefore only interpreted in terms of broad sources that produce anomalies with amplitudes much greater than 3 mGal, and surficial sources are not studied.

The margins of the batholith complex were interpreted from the horizontal gradient magnitude of the new gravity map. This method is based on the principle that for the case of vertical boundaries, the gravity field's largest magnitude gradients in the horizontal direction are located over the edges of density contrasts (Cordell, 1979; Blakely and Simpson, 1986). The method is therefore useful for detecting geologic contacts. Here, the horizontal gradient magnitude was calculated from the new gravity map after upward continuation of 5 km, in order to attenuate effects of near-surface sources (Fig. 9). The batholith complex margin, assumed to be vertical, was interpreted from the result, and the location is similar to the previously interpreted extent for the western, southern, and northern portions of the complex (Fig. 3) (Plouff and Pakiser, 1972). The interpreted eastern margin of the batholith complex is extended significantly to the east, such that it nearly reaches the western margin of the San Luis Basin (cf. Figs. 3 and 7). It is possible that interference from other anomaly sources, such as the San Juan sag, may be present in the new gravity anomaly and horizontal gradient maps. Nevertheless, the interpreted areal extent of the complex is 8200 km².

The vintage analog seismic refraction profiles (A–A′; Fig. 1) were reinterpreted to provide subsurface constraints on gravity models made from the new gravity map. First break times of refracted and reflected phases were picked from analog paper records provided by C. Prodehl (2007, written commun.), and modeled using MacRay (Luetgert, 1992), a ray-tracing forward-modeling program (Fig. 10). The picking error of phase arrival times was estimated to be 0.2–0.3 s, based on resolution of the original records. Picks were not significantly different from those originally made by Prodehl and Pakiser (1980).

Gravity and seismic refraction data were modeled in an integrated fashion along profiles A–A′, B–B′, and C–C′ (Fig 5). The fundamental assumption underlying the following discussion is that the source of the gravity low is the same as the source of the seismic low-velocity zones. Whereas not specifically addressed by previous studies, this assumption makes sense, as the source body interpreted by gravity data analysis (Plouff and Pakiser, 1972) occupies the same portion of the upper crust occupied by the sources of seismic low-velocity zones (Prodehl and Pakiser, 1980), and low seismic velocities are normally correlated with low densities. This common source is interpreted to be a silicic Oligocene batholith complex. Because seismic refraction is a useful method for detecting broad variations of velocity with depth, the seismic modeling was used to define the top and bottom of the complex, and the gravity horizontal gradient magnitude map (Fig. 9) was used to define its lateral boundaries in the models. Other assumptions used to make the models presented here are that the complex is entirely below an elevation of 2286 m (7500 ft), which the seismic model supports (see following), and that its density contrast with country rocks remains constant with depth.

The first step was to develop a crustal P-wave velocity model that satisfied the first break picks made from the seismic data along profile A–A′, a short distance to the east of the central caldera cluster (Fig. 7). This model was constructed from the shallow subsurface downward, adjusting layer thicknesses and velocities of refractors until the observed data were fit in a broad sense (Fig. 10). It was not possible to fit the arrival time of every first break pick, although this was not necessarily desirable due to the low spatial resolution and large possible time errors inherent in the data. Velocities averaging 4.1 km/s were used, extending from the surface to a depth corresponding roughly to sea level, where a very high velocity (6.2 km/s) refractor is interpreted. This may indicate the presence of Precambrian rocks under the eastern margin of the central caldera cluster.

Two low-velocity zones were modeled to be within the upper crust (0–20 km depth below sea level), based on traveltime modeling and the rapid terminations of two branches of refracted arrivals (Fig. 10). These zones were modeled using a variety of different velocities, and velocities ranging from 5.8 to 6.0 km/s were found to provide reasonable fits to the observed data. Attempts to model the low-velocity zones with values outside of this range resulted in refractor crossover distances and traveltimes that were inconsistent with the observed data. Because the velocity of rocks composing the low-velocity zones cannot be uniquely determined, there was a tradeoff between velocity and low-velocity zone thickness such that higher velocities result in a larger required thickness and lower velocities result in smaller required thickness. This tradeoff is expressed in terms of different batholith complex thicknesses for different velocities (Table 2). The best fit and therefore preferred seismic model used a velocity of 5.9 km/s for the low-velocity zones that are interpreted to define the batholith complex (Fig. 10).

The low-velocity zones are separated by a refractor with a velocity of 6.3 km/s. The shallowest low-velocity zone is apparent in reversed arrivals (e.g., shown by seismic arrivals from both north and south directions) and is therefore a robust feature of the model. The same is true of the 6.3 km/s refractor. However, the deeper low-velocity zone is only detected in seismic data from the northern shot, due to a lack of receivers south of the southern shot (Fig. 10), and so is a less robust feature, despite providing a reasonable fit to the observed data. Deeper portions of the crust and upper mantle are not well imaged by the seismic data. The gravity model along profile A–A′ was constructed to match the upper crustal configuration indicated by the seismic model (see following discussion; Fig. 11).

For profile B–B′, the first step of the gravity modeling was to define the regional gravity field, because this profile is coincident in location with the model of Roy et al. (2004), wherein the regional gravity field and interpreted source were first defined. The same source of low densities in the mantle was used here: an ∼300-km-wide, ∼200-km-thick zone with a density contrast of –33 kg/m³ with surrounding mantle rocks. However, Roy et al. (2004) only defined the regional field in an east-west direction. The upward continued gravity data were used as a rough guide to extend the interpreted regional source to the south and northwest (Fig. 5) to profile C–C′. The northern extent of the regional source (end of A–A′) was more difficult to define because the regional gravity low over the SJVF merges with a broader regional gravity low over much of the Rocky Mountains of central Colorado. The continuation of the gravity low to the north of A–A′ may be related to a number of different sources, including crustal thickening under the Rocky Mountains of central Colorado, and therefore does not necessarily indicate a continuation of low mantle densities north of A′. Given this uncertainty, the north edge of the low densities in the mantle was initially left loosely defined along A–A′.

The next step was to use the top and bottom contacts of the batholith complex interpreted by seismic modeling in the gravity model along A–A′, using the lateral position of the batholith complex edges from the horizontal gradient magnitude analysis (Fig. 9). These contacts were extended to B–B′ and C–C′ (Fig. 11). Modeling of the three profiles together allowed an overall best fit density contrast for the batholith complex to be determined, as well as variations in thickness away from A–A′. Once the gravity signature of the complex was determined across all three profiles, the location of the northern edge of the mantle low-density zone was adjusted until the model along A–A′ fit the observed data as shown in Figure 11. The entire process of using the seismic model to constrain the gravity interpretation was repeated using the range of reasonable velocities for the batholith complex (Table 2), although only the preferred models are shown here (Fig. 11). The seismic refractor that is between the two low-velocity zones was included in all the gravity profiles because it is a robust feature of the seismic model, although this is not required by the gravity data.

Neither the gravity nor seismic refraction data are capable of resolving changes of density and/or velocity with depth within the batholith complex. For the preferred model, the computed density contrast is a vertically averaged value of –80 kg/m³. This corresponds to a batholith complex density of 2620 kg/m³ if a background upper crustal density of 2700 kg/m³ is assumed (only the density contrast can be uniquely determined). This contrast can reasonably be varied from –60 to –110 kg/m³ for seismic velocities of 6.0 and 5.8 km/s, respectively (Table 2). All modeling runs showed the batholith complex to be significantly thicker under the western portion of the SJVF than under the central and eastern portions (Fig. 11). The preferred model shows a maximum thickness of ∼20 km under the western SJVF that decreases eastward to <10 km under the eastern SJVF. The model that produces a maximum thickness of 23 km, corresponding to a density contrast of –60 kg/m³, is consistent with the thickness estimate of Plouff and Pakiser (1972). However, the thickness estimate range presented here is generally thinner than the previous estimate (Plouff and Pakiser, 1972), due to the modeling approach taken here that resulted in a smaller amplitude gravity low. Assuming an area of 8200 km² (from the horizontal gradient magnitude analysis; Fig. 9) and using thickness estimates from the gravity modeling, the volume of the complex is estimated to be 82,000–130,000 km³ (Table 2).

The Southern Rocky Mountain volcanic field is generally considered to be a manifestation of the Cenozoic ignimbrite flare-up that affected much of southwestern North America. The flare-up is characterized by a high volume of erupted intermediate and silicic volcanic rocks, the latter typically associated with nested caldera collapse and regional emplacement of far-traveled ash-flow tuffs and associated deposits (Lipman, 2007, and references therein). The subarea of the SJVF, those erosional remnants constituting the San Juan Mountains, has been interpreted to preserve the predominantly volcanic record of mainly Oligocene high-power magmatic input (Lipman, 2007). This magmatism overlaps in time with the transition to extensional tectonism of the Rio Grande rift. Regardless of differing opinions of tectonic controls on mantle source materials, mechanism of melting, magma ascent, and differentiation models, the general consensus remains that Oligocene supervolcanoes of the SJVF reflect an almost unprecedented input of mafic magma into the crust, a necessary thermal requirement to generate the estimated 40,000 km³ of erupted materials (Lipman, 2007). Only the Altiplano-Puna volcanic complex in the central Andes (de Silva, 1989; Lindsay et al., 2001; Schmitt et al., 2003; de Silva et al., 2006; de Silva and Gosnold, 2007) provides a reasonable late Cenozoic to Quaternary corollary with respect to eruptive style, petrologic diversity, and scale (although total erupted volumes are less, ∼15,000 km³).

Given the extraordinary volumes of erupted magma, much of it cataclysmically erupted from multiple caldera sources, the nature of associated subvolcanic magma chambers solicits considerable discussion and debate. Constraints on volume and geometry are critical to these discussions. We focus this discussion on the geophysically observable characteristics of the subsurface environment of the SJVF, particularly crystallized magma chambers that coalesced as a large batholith complex.

Total volume estimates for the subvolcanic batholith complex of the SJVF subarea of the southern Rocky Mountain volcanic field range from 82,000 km³ to 130,000 km³, consistent with estimates of total intrusive volume of ∼100,000 km³ (Lipman, 2007). Taking into account an erupted volume of ∼40,000 km³, ratios of intrusive volumes (presumed to constitute magma chamber residuum following a major eruption) to volcanic volumes range from 2:1 to ∼3:1. These values are significantly lower than bounding estimates of 10:1 proposed for continental magmatic systems (Smith, 1979; Crisp, 1984). Lower values appear to be more reasonable for large systems where ignimbrite eruptions may drain the upper several kilometers of any given magma chamber (White et al., 2006; Lipman, 2007), as appears to be the norm in the Southern Rocky Mountain volcanic field. However, our ratios based on calculated batholith complex volumes may reflect minimums, because the lateral extent of the batholith complex is based on horizontal gradient trends of gravity data that are filtered to reflect the effect of the main mass of the batholith complex, excluding intrusions that almost certainly underlie calderas outside of the main gravity low in the northeast and southeast San Juan Mountains (Fig. 9). Perhaps more important, relatively deep and dense rocks of or related to the batholith complex may not be geophysically well imaged, meaning that they are not included in the volume estimates presented here.

Our underlying assumptions are that the source of the gravity low is the same as the source of the seismic low-velocity zones and that the density of the interpreted batholith complex does not change significantly with lateral position. Integrated modeling can thus be used to infer the batholith complex's thickness variations in addition to volume. The thickest part (15–23 km, Table 2) of the batholith complex underlies the western SJVF and the western caldera cluster, where there are four nested calderas on the surface. The complex is thinner under the central caldera cluster, with seven calderas, showing that there is not a direct relationship between batholith complex thickness and mapped calderas. We interpret this as reflecting the integrated temporal and spatial development of subvolcanic magma chambers and their crystallized plutonic equivalents over protracted periods of time (millions of years), as opposed to reflecting plutonic residuum associated with any given caldera. This is consistent with interpretations that the batholith complex was characterized by accumulation of magma after volcanism ended, and that the batholith complex remaining today is not necessarily representative of magma chambers that existed during any specific caldera-forming eruption (Lipman, 2007). A more direct relation between batholith complex thickness and caldera clustering would be expected if the complex represented fossil magma chambers that were tapped by surficial volcanism.

Despite the fact that the batholith complex is shown here to have constant physical properties and simple boundaries, it almost certainly was assembled over several million years (e.g., Lipman, 2007), is composed of rocks emplaced in multiple episodes with different densities, and has complex boundaries. The simple models presented here are only approximations of its properties and dimensions, which are chosen to be consistent with the seismic model. For example, batholiths are thought to become more mafic in composition and thus denser with increasing depth, although with no geophysical constraints on increasing density or velocity with depth, any models with this characteristic would be quantitatively arbitrary and were not attempted here. Further, the top and bottom of the geophysically imaged batholith complex may not correspond well to boundaries that would be picked based on petrologic analyses. Despite these limitations, the overall (i.e., vertically averaged) densities (2590–2640 kg/m³) determined for the batholith complex, assuming an upper crustal background density of 2700 kg/m³, are consistent with published values for granite densities (Daly et al., 1966; Telford et al., 1990), as well as density measurements on granite outcrops of the SJVF (Popenoe and Steven, 1969; Popenoe and Luedke, 1970; Grauch and Hudson, 1987).

The interpreted density range of 2590–2640 kg/m³ assumes an overall upper crustal density of 2700 kg/m³, a value that is commonly assumed for gravity studies but may or may not be valid in the SJVF region. The vertically averaged density contrast between the batholith complex and surrounding rocks, however, can be uniquely determined. The range of permissible values, –60 to –110 kg/m³, particularly the preferred value of –80 kg/m³, compares favorably to the few other studies that have examined the geophysical expressions of large batholiths and have included joint seismic and gravity modeling. The inferred batholith under the Mogollon-Datil volcanic field of New Mexico and Arizona (Fig. 1) is interpreted to have a density contrast of –100 kg/m³ (Schneider and Keller, 1994). The Boulder batholith of Montana, similar in dimensions and interpreted thickness to the SJVF batholith complex (Biehler and Bonini, 1969), was found by joint seismic reflection and gravity modeling to have a density contrast of –80 kg/m³ (Vejmelek and Smithson, 1995).

Although the main source of the gravity low is interpreted to be an Oligocene batholith complex, elsewhere in the southern Rocky Mountain region Precambrian granitic intrusions are interpreted to cause gravity lows (Quezada et al., 2004). As Precambrian granites crop out on the northern and southwestern margins of the SJVF (Fig. 2) and may floor the magmatic complex, significant contributions to the gravity low from low-density Precambrian granites buried under the SJVF cannot be ruled out, especially because their distribution at depth is not known. However, this possibility seems unlikely to us, given that Precambrian granites are normally not directly associated with gravity lows where they crop out in the SJVF region.

Another assumption in the gravity modeling and the interpreted thickness distribution is that the density of the batholith complex does not vary vertically or laterally. Although batholith density may reasonably be expected to vary vertically (Miller and Miller, 2002), it is difficult to model these variations in a meaningful way because little is known of the batholith complex composition at depth. Lateral density variations are also possible, as observed dramatically within the Sierra Nevada batholith in California (Oliver, 1977; Oliver et al., 1993). Such variations within the SJVF batholith complex would change the interpretation of lateral thickness patterns. However, this is difficult to address quantitatively because the batholith complex cannot be sampled for densities, and assignment of variable densities would be arbitrary and unsubstantiated. Existing density measurements on Oligocene pluton outcrops are confined to the western portion of the study area and further may not be representative of the deeper batholith complex (Lipman, 2007).

The batholith complex underlying the SJVF is expressed seismically by relatively low velocities compared to surrounding upper crustal rocks. Numerous studies have interpreted seismic low-velocity zones to indicate magma chambers under active caldera systems of similar scale to the SJVF (de Silva, 1989; Chmielowski et al., 1999; Zandt et al., 2003; de Silva et al., 2006; Bachmann et al., 2007), although this study is one of few to interpret low velocities for a batholith complex (i.e., a crystallized magma chamber or chambers). A similar interpretation was made for the batholith under the Mogollon-Datil volcanic field (Schneider and Keller, 1994) that is comparable in age to the SJVF, and also along the western margin of the Rio Grande rift (Fig. 1). A 1–2-km-thick refractor (e.g., high-velocity zone) between two low-velocity zones was found to be a robust feature of the seismic model along A–A′. The geologic significance of the refractor is unknown, although it has a higher velocity than low-velocity zones that are above and below. Its presence and geometry are not constrained by the gravity data and modeling, although gravity models along B–B′ and C–C′ show it in order to maintain consistency with the models along A–A′ (Fig. 11). Although it is compelling to try to interpret this high-velocity zone within the broader low-velocity mass as a manifestation of magmatic processes (e.g., resulting from a two-stage process of magma emplacement), any such interpretation would be poorly constrained and largely conjecture.

The timing of batholith complex emplacement is not directly addressed by the geophysical analyses presented here. However, the relatively old calderas of the northeastern and southeastern SJVF (Fig. 2) are outside of the area of interpreted batholith extent (Fig. 9). While speculative, we suggest that batholith complex emplacement was more closely associated with the formation of the younger central and western caldera clusters.

A relationship may exist between batholith complex thickness and the elevation of the overlying ground surface. The thickest part of the complex is interpreted to underlie the western SJVF, spatially coincident with the highest elevations of the San Juan Mountains (Figs. 4 and 11). This observation is consistent with a previous hypothesis that thick, low-density batholiths may cause uplift of overlying rocks on the order of hundreds to several hundreds of meters (Lipman, 1988, 2007). However, we consider this to be speculative, and do not address it further.

A new approach to gravity data processing, one that used geologically appropriate densities and digital topography to represent areas of high elevations, was used to compute a new complete Bouguer anomaly–like gravity map of the San Juan Mountains region. Traditional methods of gravity data processing cause interpretation problems there due to high elevations being composed of rocks of the SJVF that have densities significantly lower than the standard Bouguer correction density (2670 kg/m³). A large-amplitude gravity low over the region has long been interpreted to reflect a silicic Oligocene batholith complex that is related to the SJVF. Although the processing and interpretation scheme used here resulted in a residual gravity anomaly that is 10–15 mGal smaller than that observed on a standard Bouguer anomaly map, a large gravity low remained.

The enhanced gravity expression of the interpreted batholith complex was modeled quantitatively in conjunction with reinterpreted seismic refraction data that image thick low-velocity zones in the upper crust beneath the SJVF. Assuming that the source of both the gravity low and the source of the low seismic velocities is an Oligocene batholith, a series of two-dimensional modeling runs indicated that the batholith has an overall vertically averaged density contrast of –60 to –110 kg/m³ with surrounding upper crustal rocks, an overall density of 2590–2640 kg/m³, and an estimated volume of 82,000–130,000 km³. Comparison to erupted volumes yields ratios of intruded to extrusive magma volumes of ∼2:1 or 3:1, although these likely represent a minimum because the deepest and most dense parts of the batholith complex may not be geophysically well imaged. The interpreted density contrasts, densities, and inferred seismic expression of the batholith are consistent with other similar intrusions in the western United States.

The batholith complex is thickest under the western SJVF, as much as ∼20 km in the best fit model, and thins to <10 km under the eastern SJVF. The thickness pattern is not correlated with clustering of calderas on the surface, in that the largest caldera cluster does not overlie the thickest part of the complex. This is consistent with the latest petrologic interpretations that the complex continued to be modified following cessation of volcanism, and therefore is not necessarily representative of synvolcanic magma chambers.

This work was supported by the University of Oklahoma and the University of Texas at El Paso. Numerous discussions with Peter Lipman and Tom Steven helped to sharpen our arguments. Kevin Crain and Steve Holloway provided critical support to the three-dimensional gravity modeling effort, and Jefferson Chang assisted with the seismic modeling. Robbie Gries contributed gravity and density data that improved the quality of the models presented here. Reviews from V.J.S. Grauch, Rick Saltus, Chris Fridrich, Michael Cheadle, Scott Paterson, Peter Styles, and two anonymous reviewers greatly improved the quality of the manuscript. Any use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the U.S. Government.