Hydrologic models of modern and fossil geothermal systems in the Great Basin: Genetic implications for epithermal Au-Ag and Carlin-type gold deposits

The Great Basin physiographic province of the western United States is the world's second leading producer of gold and is also known for its geothermal energy production (Williams, 2002; Hofstra and Wallace, 2006). Most of the gold in this region is produced from Carlin-type deposits, pluton-related porphyry to distal disseminated deposits, and epithermal deposits (John et al., 2003; Wallace et al., 2004). This paper focuses on quantitative analysis of hydrothermal systems in north-central Nevada that generated gold deposits during three periods in different geotectonic settings: Eocene (42–36 Ma) Carlin-type Au deposits associated with calc-alkaline magmatism and extension (Hofstra and Cline, 2000; Cline et al., 2005), Miocene (17–12 Ma) low sulfidation Au-Ag deposits associated with tholeiitic bimodal magmatism and extension (John et al., 2003; John and Wallace, 2000; John, 2001), and late Miocene–Holocene (7–0 Ma) low sulfidation Au-Ag deposits and active geothermal systems associated with range-front normal faults that accommodate differential translation along the Walker Lane in southwestern Nevada (Coolbaugh et al., 2005; Faulds et al., 2005). Although inputs of heat and/or fluids of magmatic or metamorphic derivation are evident in some deposits, the genetic models for these systems call upon convection of meteoric water through underlying sedimentary, metamorphic, or igneous rocks with upflow of hot fluids along permeable high-angle faults and deposition of silica and/or Au at shallow levels. The common occurrence of gold in these systems may reflect the ability of underlying reduced carbonaceous and pyritic sedimentary rocks to shift the redox and sulfidation state of hot fluids to conditions that facilitate gold transport by sulfide complexes (Hofstra and Emsbo, 2007). Hence, an important step to increase our understanding of these systems is to quantify hydrothermal circulation in the context of the geologic framework and physicochemical constraints.

We embarked on a program of cooperative research between geologists at the U.S. Geological Survey and hydrologists at Indiana University to enhance understanding of deep groundwater flow in these systems. This work is part of a 5 yr U.S. Geological Survey project on the metallogeny of the Great Basin that seeks to understand the interplay between crustal evolution, fluid flow, and ore formation (Hofstra and Wallace, 2006). In this investigation we employed a “present is the key to the past” approach that uses information and simulations of modern aquifers and geothermal systems to constrain and inform our simulations of fossil hydrothermal systems that, respectively, produced low sulfidation Au-Ag deposits and Carlin-type gold deposits. Figure 1 and Table 1 show the locations and attributes of the systems selected for this study: Dixie Valley, Beowawe, Mule Canyon, and the northern Carlin trend. During the course of this investigation a number of regional and system scale geologic cross sections and three-dimensional algorithms were modified and the role of high-permeability fault zones was evaluated (Gao et al., 2003; Person et al., 2005), and a new parallel multi-component, finite element model was developed to simulate crustal-scale fluid flow, heat, isotope, and silica transport (Person et al., 2007). The new model was used to test several hypotheses, such as the effects of magmatic intrusion on ore genesis (Banerjee et al., 2007), that have enhanced our understanding of these systems. The purpose of this paper is to present an overview of our first steps to quantify deep groundwater convection in these systems using cross-sectional models. Our ongoing work is focused on evaluating the impact of inputs of heat and/or fluids from deep sources on these systems in a three-dimensional framework using a new, parallel three-dimensional hydrothermal flow model we have recently developed (Wang et al., 2007).

We first review shallow heat flow data and temperature logs from geothermal wells completed within the Beowawe and Dixie Valley geothermal fields. A regional-scale (100 km) hydrogeologic model is then presented to evaluate the possibility that topographically driven regional groundwater flow systems within the Paleozoic carbonate aquifer system could explain the occurrence of the Battle Mountain heat flow high and the Eureka heat flow low (Blackwell, 1983; Fig. 2). Next, we present a more local model of the present-day Beowawe geothermal field, which incorporates a series of subvertical, highly permeable faults connected to each other through a thin (~40 m), highly permeable aquifer within the Paleozoic carbonate unit. Using a similar range of permeability consistent with field conditions observed today at Beowawe, we present a cross-sectional model of the Miocene Mule Canyon low sulfidation epithermal Au-Ag system and the Eocene Carlin-type Au system in the northern Carlin trend. Our quantitative reconstructions of these fossil flow systems are constrained by comparisons between computed and observed temperatures, stable isotopic data, silica concentrations, apatite fission track ages, and total volume of silica precipitated at the land surface.

The hydrogeology of the Great Basin has been studied in detail because of municipal and agricultural water demands, mine dewatering, and nuclear waste isolation issues (Harrill and Prudic, 1998). Most agricultural related hydrogeological studies have focused on characterizing shallow flow conditions (<500 m) within alluvial aquifers (e.g., Prudic and Herman, 1996). Municipal water supply, mine dewatering, and nuclear waste isolation studies have considered groundwater flow conditions within the deeper (>2 km) Paleozoic carbonate aquifer system (Welch and Bright, 2007; Maurer et al., 1996; Belcher, 2004). Plume and Carlton (1988) indicated that important aquifers of the Great Basin include shallow alluvial aquifers of Quaternary–Tertiary age (to 400 m thick), and Paleozoic carbonate rock aquifers (to 5 km thick). Paleozoic siltstones, shales, and Mesozoic–Cenozoic volcanic units generally are characterized as lower permeability confining units, although some welded tuffs have significant fracture permeability and function as locally important aquifers (Blankennagel and Weir, 1973; Winograd and Thordarson, 1975). The Paleozoic carbonate aquifer is composed of upper and lower permeable members separated by a Mississippian–Devonian confining unit (Plume and Carlton, 1988; Dettinger et al., 1995). Aquifer tests within carbonate rocks (Plume and Carlton, 1988; Dettinger et al., 1995) yield hydraulic conductivity values between 0.03–10 m/day (10^−11.1–10^−9.1 m²). This is comparable to values reported by Winograd and Doty (1980) in southwestern Nevada in the vicinity of the Nevada Test Site. The bulk hydraulic conductivity of carbonate rocks in the area of Post-Betze Mine from the dewatering of the mine is ~30 m/day (Hydrologic Consultants, 1999). Permeability within the Paleozoic carbonates is controlled by fracturing, faulting, and solution features associated with both Cenozoic and Paleozoic episodes of karstification, whereas the high permeability of the siliciclastic aquifers is mainly due to fractures (Maurer et al., 1996). Siliciclastic confining units are reported to have a hydraulic conductivity of between 0.0003–0.1 m/day (10^−15.5–10^−13.0 m²). The hydraulic conductivities of Quaternary pluvial, fluvial and/or alluvial, and colluvium units range between 0.15 and 610 m/day (10^−10.4–10^−6.8 m²). Tertiary basin fill is reported to have three orders of magnitude lower hydraulic conductivities (10^−9.5 m²).

Movement of groundwater in the Great Basin is influenced by a combination of topography, climate, and geology. Groundwater is believed to move through permeable zones under the influence of hydraulic gradients from areas of recharge to areas of discharge (referred to herein as topography-driven flow; Toth, 1963). In alluvial basins, water-table patterns and flow directions mimic topography (e.g., Prudic and Herman, 1996). Groundwater recharge to the shallow alluvial aquifers occurs in topographically elevated alluvial fans fed by surface water runoff derived from losing streams issuing from adjacent ranges especially during periods of spring snow melt. Discharge occurs in lowlands as base flow to streams and playas. Deep regional groundwater flow through consolidated rock aquifers is unconstrained by local topographic or drainage features and is driven by hydraulic gradients that are continuous over tens to hundreds of kilometers. The principal consolidated-rock aquifer is a thick sequence of Paleozoic carbonate rock that extends throughout the subsurface over much of eastern and central Nevada and forms regional multibasin flow systems (Dettinger et al., 1995; Harrill and Prudic, 1998). There are five regional flow cells within the carbonate aquifer system in Nevada: the Upper Humboldt River, Railroad Valley, Death Valley, Bonneville region, and Colorado River region (Fig. 2). The solid blue lines denote the boundary of these flow systems. The black contours denote the freshwater equipotential contours, and the black dashed flow lines with arrows depict the inferred groundwater flow directions (Fig. 2) (assuming topography-driven flow). Recharge to the regional groundwater flow systems primarily originates in mountains; the magnitude of recharge is strongly controlled by elevation within mountain blocks (Maxey and Eakin, 1949). Recharge estimates for Beowawe and Dixie Valley range between 5.5 and 0.6 cm/yr in uplands with no recharge occurring in most valleys (Harrill and Hines, 1995; Olmsted and Rush, 1987). Discharge from the regional groundwater flow systems is manifested by large springs and, in some areas, extensive wetlands (Thomas et al., 1986). Most shallow groundwater discharge is consumed by evapotranspiration. Regional springs are characterized by (1) high discharge rates, >~550 m³/day, at relatively constant volumes (Harrill et al., 1988); (2) elevated temperatures, >27 °C, assumed to be indicative of deep circulation (Mifflin, 1968); and (3) by water-chemistry indicative of flow through consolidated carbonate rock (Mifflin, 1968; Steinkampf and Werrell, 2001). Regional groundwater flow in the deeper carbonate aquifer system results in significant interbasin fluid transfer (Harrill and Prudic, 1998; Nichols, 2000; Welch and Bright, 2007). These regional flow systems are responsible for significant lateral transport of heat, but probably do not account for the occurrence of modern geothermal features such as Dixie Valley or Beowawe, which appear to require the addition of magmatic heat or fluid circulation along deep-seated faults (Coolbaugh et al., 2005). The Great Basin is known for two broad heat flow anomalies (Lachenbruch and Sass, 1978; Fig. 2): the Battle Mountain heat flow high (heat flow >104 mW/m²) and the Eureka heat flow low (heat flow <63 mW/m²; Sass et al., 1971). These broad heat flow anomalies are overprinted with dozens of more localized geothermal fields (Benoit and Butler, 1983; Hulen et al., 1990) that have measured shallow heat flow anomalies as high as 4000 mW/m² (Olmsted and Rush, 1987). Blackwell (1983) called on crustal thinning and associated enhanced basal heat flow to explain the Battle Mountain heat flow high. This study, however, did not consider convective effects associated with groundwater flow. Figure 2 shows that regional high heat flow anomalies (Battle Mountain) and Paleozoic carbonate aquifer system discharge areas are not well correlated. It is important to note that the Battle Mountain heat flow high is just to the west of the termination of the Paleozoic carbonate aquifer system (Fig. 2) and thus the deep aquifer system cannot directly contribute heat to this geothermal feature. However, the Eureka heat flow low is situated near regional flow systems divides in recharge areas of the Death Valley, Railroad Valley, and Colorado River flow systems. We evaluate whether deep groundwater flow systems are responsible for regional heat flow anomalies later in this paper using numerical models.

While about half of the Basin and Range geothermal systems, such as Coso, are clearly influenced by nearby magmatic activity (Monastero et al., 2002; Person et al., 2006), many are not (Coolbaugh et al., 2005). Of the nine high-temperature (>200 °C) geothermal systems in the northern Basin and Range, five have no evidence of recent magmatic activity (Benoit and Butler, 1983). The regional transtensional setting appears to provide the fracturing and permeability needed at sufficient depths to form such high-temperature geothermal systems without magmatic heat (Coolbaugh et al., 2005). Inspection of heat flow maps presented by Coolbaugh et al. (2005) indicates that zones of high shallow temperature gradients (>100 °C/km) only have a lateral footprint of 1–3 km. However, a comparison of the volumetric fluxes from thermal springs (in excess of 35 °C) and consumptive use from phreatophytes suggests that geothermal discharge is a small portion of the shallow groundwater budget across the Basin and Range Province (D. Prudic, 2008, personal commun.).

Generic (López and Smith, 1995) and site specific (Fairley and Hinds, 2004; McKenna and Blackwell, 2004) numerical models demonstrate that deep, permeable fault zones are needed to account for such small heat flow anomalies. Using transient hydraulic fluctuations associated with Earth tides combined with numerical modeling, Bredehoeft (1997) estimated the permeability of subvertical faults cutting the carbonate aquifer system in the Basin and Range near Yucca Mountain to be ~10⁻¹⁰ m². There is geothermal and isotopic evidence (i.e., elevated temperatures along faults as well as anomalously young groundwater ages) that faults in this region act as conduits connecting the deeper carbonate aquifer with overlying hydrogeologic units (Winograd and Pearson, 1976; Bredehoeft, 1997; Rose and Davisson, 2003). An important question we consider herein is whether groundwater circulation in these geothermal systems is generated by topography-driven lateral flow through permeable aquifers with discharge along range-front faults, as suggested by McKenna and Blackwell (2004), or whether it is induced by natural convection cells within the plane of range-bounding fault zones (López and Smith, 1995).

We begin our study with a field examination and review of the literature on Dixie Valley because it is one of the largest geothermal power producers in the Great Basin and has been the site of numerous studies. We then focus on the Beowawe geothermal system and the Miocene hydrothermal system that formed the Mule Canyon low sulfidation Au-Ag deposit because they have been significant producers of energy and gold, respectively (Table 1), they have been studied in some detail (e.g., John et al., 2003, and references therein) and they are located close to one another in Whirlwind Valley (Fig. 1). Their close proximity and similar settings made them ideally suited for our approach to use the constraints on modern systems to inform models for fossil systems. The sedimentary rock units in the subsurface of the Northern Nevada Rift also are similar to those in the Carlin trend, which is the largest producer of gold in the Great Basin and the site of numerous investigations. We then use the results of our work on the younger systems to inform development of models for the Eocene hydrothermal system in northern Carlin trend. Figure 1 shows the locations of these systems in northern Nevada, and Table 1 provides a synopsis of their attributes and parameters that constrain fluid flow models.

Heat flow geophysicists, hydrogeologists, and economic geologists have come to appreciate the important association between high heat flow and hydrothermal activity along tectonically active Basin and Range fault systems that accommodate differential translation along the Walker Lane (Faulds et al., 2004, 2005; Coolbaugh et al., 2005). This geotectonic setting has existed for as much as 7 m.y. and has been called upon to explain the spatial association of young low sulfidation Au-Ag deposits and active geothermal systems with range-front faults. An example relevant to this study is the Pleistocene Dixie Comstock Au-Ag deposit on the eastern margin of the Stillwater Range (Vikre, 1994) and the active Dixie Valley geothermal system (Blackwell et al., 2002) located along the same range-front fault. Other prominent examples include the 4–1 Ma Crofoot-Lewis Au-Ag deposit on the western margin of the Kamma Mountains near Winnemucca, Nevada (Ebert and Rye, 1997), and the 5–2 Ma Florida Canyon Au-Ag deposit and active Humboldt House geothermal system on the western margin of the Humboldt Range (Samal et al., 2005; Samal and Fifarek, 2007; R. Fifarek, 2005, personal commun.). The Florida Canyon Au-Ag deposit is notable in that it is hosted in sedimentary rocks and has some features in common with Eocene Carlin-type gold deposits (Hofstra et al., 2005), such as disseminated Au-bearing pyrite and pervasively silicified limestone (jasperoid). These examples suggest that an improved understanding of modern geothermal systems can provide insights into older hydrothermal systems that generated significant gold deposits.

The Dixie Valley geothermal field is located between the Stillwater and Clan Alpine Ranges. Hydrothermal activity occurs along a deeply penetrating (at least 3 km deep) range-bounding fault (Okaya and Thompson, 1985; Fig. 3) that crops out along the base of the Stillwater Range on the western margin of Dixie Valley. It hosts a 66 MW geothermal power plant, one of the largest in the Great Basin. Precipitation varies within adjacent ranges between ~12 and 33 cm/yr (Harrill and Hines, 1995). Recharge in the ranges varies between 0 cm/yr below an elevation of 1524 m to a maximum of 5 cm/yr above 2740 m (Harrill and Hines, 1995). Geophysical surveys and geologic mapping have defined the structural features and thicknesses of basin fill associated with this geothermal field (Allis et al., 1999; Smith, 1968; R. Smith et al., 2002). Aquifer tests by Morin et al. (1998) report transmissivity within the range-bounding fault plane to range between 10^−3.4–1.2 m²/min (10^−14.1–10^−9.7 m² assuming a fault plane width of 10 m). These aquifer tests produced a damped oscillatory pressure response indicative of high-permeability conditions (Shapiro, 1998). Hickman et al. (1998), Barton et al. (1998), and Blewitt et al. (2002) argued that the high-permeability conditions along this range-bounding fault are due to the stress state and strain rates of the area. Temperatures reach 279 °C at a depth of 3 km in some of the hottest geothermal wells (Blackwell et al., 2002; Fig. 3A). McKenna and Blackwell's (2004) hydrothermal model was able to reproduce observed temperatures within the Dixie Valley field assuming topography-driven recharge within the adjacent Stillwater Range and discharge occurring along the range-bounding fault (Fig. 3C). No natural convection cells are apparent in these models, although there must be some fluid density effects on flow. A problem with this model is the lack of modern precipitation in this arid region and the Pleistocene age of the groundwater in the system (Nimz et al., 1999; Janik, 2002). McKenna and Blackwell (2004) did not find evidence of any Holocene water, suggesting that shallow recharge does not contribute to this geothermal system. One would expect a mixture of both Pleistocene and Holocene water if the Dixie Valley system were driven by topography-driven flow. Our analysis of geothermal exploration well temperature profiles indicates that there is a region of shallow low heat flow associated with geothermal exploration wells 52–14, 66–16, and 72–13 (Fig. 3A) near the southern terminus of Dixie Valley that has been overlooked. The low heat flow in these wells suggests that recharge into the range-bounding fault is occurring within the basin lowlands. M. Coolbaugh (2007, personal commun.) noted additional wells with low heat flow. The presence of both a heat flow low and a heat flow high along the range-bounding fault is suggestive of a natural convection groundwater flow system in the plane of the fault. This example helped us realize the potential importance of free convection in high-permeability faults and prompted us to modify our gridding algorithms and numerical models to represent this. The relative importance of groundwater flow driven by water-table gradients and natural convection are discussed later.

The Beowawe geyser field is in a surface hydrologic drainage basin that encompasses much of Whirlwind Valley, the Argenta Rim, and the Malpais Range (Zoback, 1979; Fig. 4A). It is situated within the Battle Mountain heat flow high and is hosted in Cenozoic volcanic rocks that are presumably underlain by the Paleozoic rocks of the Upper Humboldt River regional flow systems. Prior to the installation of a 16 MW power plant, Beowawe hosted more geysers than any other geothermal field in the United States outside of Yellowstone National Park. Permeability measurements along the Malpais fault vary between 1000 and 200 Darcy-ft (10^−10.6–10^−11.3 m² assuming a fault width of 10 m; Faulder et al., 1997). Although it contains altered and mineralized rocks with anomalous concentrations of Au and Ag, significant mineral resources have not been identified (Cole and Ravinsky, 1984; Struhsacker, 1986). The Beowawe geothermal field is situated along the east-northeast–striking Malpais fault, which cuts across the north-northwest–trending Northern Nevada Rift (Fig. 1 in John et al., 2003). Faulder et al. (1997) described complex, fault-controlled, three-dimensional flow system patterns here that may be similar to the situation at Dixie Valley (described above). Groundwater discharge within the geothermal field is confined to the Malpais fault zone, where as much as 10¹¹ kg of siliceous sinter has been deposited at the surface (Rimstidt and Cole, 1983). The field is located along the contact between consolidated rock and basin fill ~50 m above the floor in Whirlwind Valley. One reason for the occurrence of the geyser field may be that the basin fill contains a thick unit of fine-grained deposits (lacustrine) that would inhibit groundwater flow across the basin fill. Another reason may be deep circulation of water along the southwest-northeast–trending fault that coincides with the field (David Prudic, 2008, personal commun.; Fig. 4A). Prior to development, discharge from the Beowawe geyser field was ~61 L/s (Waring, 1965). While the Beowawe geothermal field is characterized by remarkably high temperatures (215 °C at a depth of only 2.9 km; Zoback, 1979) and shallow heat flow (>2000 mW/m²; Olmsted and Rush, 1983), it is apparently not magma related, based on the absence of Quaternary volcanic rocks (Coolbaugh et al., 2005) and lack of enrichment in mantle-derived ³He (White, 1998). In addition, stable isotopic analyses of modern geothermal waters from Beowawe have a meteoric composition indicating little fluid-rock isotopic exchange (John et al., 2003). Tracer tests have been used to determine the connectivity of different wells within the geothermal field (Rose et al., 1995). Zoback (1979) used dissolved silica concentrations, geyser discharge rates, and the volume of sinter deposits to estimate that the sinter deposits were formed in ~200 k.y. Shallow (Olmsted and Rush, 1987) and deep (Smith, 1983) geothermal wells in the vicinity of the Beowawe geothermal field (Fig. 4B) display a concave temperature profile. These data taken together are consistent with a fault-controlled hydrothermal system (López and Smith, 1995).

While local hydrologic budget calculations led Olmsted and Rush (1987) to conclude that hydrothermal circulation associated with Beowawe is confined to the local surface drainage of Whirlwind Valley, their surface heat flow survey suggests that much of this watershed is underlain by an upflow zone (i.e., heat flow >80 mW/m²). The recharge area to this geothermal system should be characterized by a heat flow low (<60 mw/m²; e.g., Smith and Chapman, 1983), which appears to be absent in Whirlwind Valley. We propose that recharge to the Beowawe geyser system comes from infiltration along permeable faults near Mule Canyon gold deposits on the west side of Whirlwind Valley. Our cross-sectional model evaluates this scenario. We further hypothesize that these recharge zones are not readily apparent because they are relatively short-lived (<10 k.y.) and thus do not create a broad lateral thermal footprint and have remained undetected.

Middle Miocene low sulfidation epithermal Au-Ag deposits in the northern Great Basin are associated with a period of bimodal basalt-rhyolite magmatism and east-northeast extension related to the Yellowstone mantle plume (Wallace and John, 1998; John, 2001; Ponce and Glen, 2002). Regional heat flow was probably high. Anhydrous, reduced, bimodal basalt-rhyolite magmas erupted along several north-trending normal faults, such as the Northern Nevada Rift, which is manifest by a 500-km-long aeromagnetic anomaly that extends from near the Oregon-Nevada border to southeastern Nevada (Zoback et al., 1994). Its surface trace is marked by middle Miocene volcanic centers and epithermal Au-Ag and Hg deposits. Mule Canyon is one of the few significant Au-Ag deposits of the group associated with a mafic volcanic center; most of the others are associated with rhyolite domes (John et al., 2003).

The Mule Canyon Au-Ag deposit contains ~1 Moz (~31 t) of Au and is located ~8 km northwest of the Beowawe geothermal field near the west edge of the Northern Nevada Rift (Fig. 1). The north-northwest–striking Muleshoe normal fault was active during formation of the rift and served to localize a basalt-andesite eruptive center composed of dikes and flows dated at 16.4–15.8 Ma and the Au-Ag deposit dated at 15.6 Ma (John et al., 2003). Prior to the eruption of volcanic rocks, the area accumulated lacustrine sediments. Some of the older flows are interbedded with tuffaceous sediments. The presence of palagonized flows and hydrated glassy flows with vague pillow structures suggests eruption into water. Hydrothermal activity was ~200 k.y. after the last volcanic eruption at Mule Canyon, such that the dikes and flows that host the deposit were cold by the time the hydrothermal system was active; eruption of mafic magmas farther to the north may have supplied heat for the hydrothermal system (John et al., 2003). However, geophysical data suggest that larger magma bodies were present at mid-crustal depths (Ponce and Glen, 2008) that may have still been hot. The evidence presented by John et al. (2003) also places important constraints on the pressure-temperature conditions and nature of the fluids in the system (Table 1). The deposit formed at shallow depths (<150 m) from dilute, near neutral pH, reduced, H₂S-rich, CO₂-bearing water. Fluid inclusions record peak temperatures of 230–265 °C at lithostatic pressures, but most of the silica in the deposits formed at temperatures <200 °C under hydrostatic conditions during periods of boiling. While silica sinter may have formed at the paleosurface, none is preserved. Stable isotope data indicate that ore fluids consisted dominantly of meteoric water that evolved by deep circulation through Paleozoic sedimentary rocks at low water:rock ratios (~1) and temperatures >200 °C. As a consequence, the δ¹⁸O values of ore fluids were shifted ~10‰ from the meteoric water line (John et al., 2003). While some of the H₂S and CO₂ in ore fluids may have come from degassing magmas at depth, the isotopic data show that most of the H₂S and CO₂ in the fluids (and probably ore metals and trace elements such as Se) were derived from marine carbonate and sulfide minerals in underlying Paleozoic sedimentary rocks (John et al., 2003).

Most of the gold produced in the Great Basin has come from Carlin-type gold deposits. Hofstra and Cline (2000), Cline et al. (2005), and Emsbo et al. (2006) provided comprehensive overviews on the geotectonic setting, attributes, and physicochemical constraints on the hydrothermal systems that formed Carlin-type gold deposits. Interpretations that frame model development are summarized in the following. These deposits are concentrated in northwest-, north-, and northeast-striking mineral belts called trends, with long and complex histories of deformation, magmatism, and hydrothermal activity (Fig. 1). Many deposits are exposed in erosional windows through antiforms in the Roberts Mountains allochthon. The Roberts Mountains allochthon is a regionally extensive Early Mississippian fold-and-thrust belt composed of a eugeoclinal Cambrian–Devonian siliciclastic slope and rise continental margin assemblage that was thrust eastward over a miogeoclinal slope and shelf carbonate assemblage. The trends and antiforms that localize the deposits are currently thought to be underlain by fault zones inherited from Neoproterozoic rifting that were reactivated in the early Paleozoic as normal faults and subsequently inverted during Antler and/or younger contractional orogenies. The deposits are broadly coeval with southward-sweeping patterns of calc-alkaline magmatism inferred to have been produced by rollback of a formerly flat subduction zone. They formed during a pulse of Eocene (42–36 Ma) northwest to west extension that dilated favorably oriented faults and channeled ore fluids upward into structural culminations.

Genetic models for Carlin-type deposits hinge on the geologic features mapped in each trend or mining district and the assortment of pale-otemperature and isotopic data gathered that constrain the source(s) of ore fluid components and corresponding fluid pathways. The minimum requirement, which is evident in every district, appears to be convection of meteoric water through reduced sedimentary rocks containing pyrite and carbon with focused discharge of hot fluids up dilatant faults and outward into permeable carbonate rocks below the Roberts Mountains allochthon. Discharge to the surface is prevented by fine-grained siliciclastic rocks that overlie the carbonates and acted as a confining unit. Gold-bearing pyrite precipitated as H₂S sulfidized iron-bearing minerals in the host rocks. In some districts there is evidence for inputs of magmatic or metamorphic fluids from below, and/or the presence of Devonian sedimentary exhalative (sedex) Au deposits, or enrichments, in the convection cells that may have been a source of Au. Isotopic evidence of fluid interaction with underlying Neoproterozoic siliciclastics is evident in one district (Getchell), likely occurred in the deeper part of the convection cells in each district, and may have been another source of Au (cf. Ilchik and Barton, 1997). The wide range of metal endowments present in each district may be a reflection of such differences.

The Carlin trend (Fig. 1) is the world's third largest gold camp and contains about half of the gold in Carlin-type deposits (100 Moz [3100 t]). The clusters of deposits in the northern, central, and southern parts of the Carlin trend (Fig. 1) generally are thought to represent separate hydrothermal systems or separate discharge areas within a larger system. The largest gold deposit, Betze-Post (40 Moz [1240 t]), the largest high-grade deposit, Meikle (7 Moz [217 t] at 24.5 g/t; Emsbo et al., 2003), and the original Carlin deposit occur in the northern Carlin trend, which was the focus of our study. The age of mineralization is constrained to ca. 39 Ma by a 39.8 ± 0.6 Rb-Sr date on late ore stage galkhaite, crosscutting relationships between mineralization and igneous dikes with ⁴⁰Ar/³⁹Ar dates that range from 40.1 to 37.3 Ma, and a pooled apatite fission track date of 37.3 Ma (Cline et al., 2005; Table 1). The age, composition, and distribution of Eocene dikes and aeromagnetic anomalies suggest that the Carlin trend is underlain by several plutons (Ressel and Henry, 2006). The largest magnetic anomaly is between the northern and central deposit clusters near exposures of the Cretaceous Richmond stock and Eocene dikes. Hickey et al. (2003a, 2003b) and Cline et al. (2005) showed that the northern end of this magnetic anomaly is approximately cospatial with a broad Eocene apatite fission track (AFT) anomaly. The Eocene AFT anomaly is superimposed on the regional pattern of Cretaceous exhumation and cooling (130–60 Ma AFT ages) obtained from areas to the north and south of the Carlin trend (Fig. 5B). Eocene reconstructions and modeling of AFT data constrain the depth of formation of the deposits to <2–3 km. The sample locations for the AFT data presented in Cline et al. (2005) are shown in Figure 1 relative to the deposits and our cross section. The spatial variability in fission track ages is due to hydrothermal activity and subsequent cooling and uplift and was used to constrain our flow models.

The Eocene landscape in the vicinity of the northern Carlin trend was characterized by moderate relief with hills and extensional half grabens filled with fluvial, lacustrine, and volcaniclastic sediments (Haynes, 2003). The northern Carlin trend was a hill on the western side of the Elko basin, which contains as much as 1500 m of fill. Mapping and geo-chronologic studies in the western part of this basin show that there was a pulse of extensional faulting and magmatism between 39.7 and 39.4 Ma that was more or less contemporaneous with ore formation. Normal faults in this basin may have been recharge areas for the northern Carlin trend hydrothermal system.

Most of the gold deposits extend over vertical intervals of <100–500 m and are localized in permeable strata within 100 m of the base of the Roberts Mountains allochthon (Fig. 5A). Such strata include karst and erosion surfaces in the carbonate platform margin, debris flows and turbidites shed from the carbonate platform into adjacent slope and basin environments, and low-angle faults and fault breccias produced by subsequent contractional and extensional orogenies (Cline et al., 2005). We therefore include a thin, high-permeability carbonate aquifer below the Roberts Mountains allochthon in our hydrogeologic model for the northern Carlin trend.

Cooling is evidenced by the dissolution of carbonate and precipitation of large amounts of quartz (or silica) to form jasperoids. Fluid inclusions in ore stage quartz from the deposits mostly yield homogenization temperatures of 160–260 °C and salinities of <3–6 equivalent wt% NaCl. They also contain a few mole percent CO₂, minor amounts of other gases (CH₄, N₂, Ar, and He) and sufficient H₂S (>10⁻² m) to transport significant concentrations of Au.

There is little isotopic evidence for magmatic water, CO₂, H₂S, or Pb in the Carlin trend. Isotopic data from the large high-grade Meikle gold deposit suggest that inputs of magmatic fluids are not required (Emsbo et al., 2003). It formed from moderately exchanged meteoric water with CO₂, H₂S, and Pb from sedimentary sources. While the northern Carlin trend has a few ore zones like North Betze and Screamer with near 0‰ pyrite that could be magmatic, most deposits have higher δ³⁴S values consistent with sedimentary sources. Likewise, although some of the hydrogen isotope data on clay minerals from Deep Star may be evidence for inputs of magmatic fluid (Heitt et al., 2003), most of the data on ore stage kaolinite (and dickite) yield calculated δD (−150‰ to −115‰) and δ¹⁸O (−7‰ to 7‰) values for water that are shifted ~10‰–20‰ from the meteoric water line (Table 1) corresponding to water:rock ratios <0.1. So, while it is possible that there were minor or local inputs of magmatic water or H₂S into the system, a significant proportion of the gold in the northern Carlin trend deposits was transported by meteoric water that evolved to become an ore fluid by reactions with reduced marine sedimentary rocks.

The geotectonic settings and attributes of these systems suggest that convection of meteoric water through reduced sedimentary rocks in response to extensional faulting, without and with magmatism, can form gold deposits. There appear to be two end-member types of such systems in the Great Basin. At one end are geothermal systems like Dixie Valley, where convection is mainly in the plane of a range-front normal fault, and Beowawe, where it is mainly along normal faults and through the high-permeability carbonate aquifer below the Roberts Mountains allochthon in the deeper part of the system. They have the smallest oxygen isotope shifts and contain small amounts of Au. At the other end is the Eocene northern Carlin trend system, where there was deep convection through carbonate and siliciclastic rocks with upward flow along faults and outward flow along the high-permeability carbonate aquifer below the Roberts Mountains allochthon at the top of the system. It has the largest oxygen isotope shifts and world-class endowments of gold. The Miocene Mule Canyon low sulfidation Au-Ag system is intermediate in character. While it had a flow path similar to Beowawe, the oxygen isotope shifts were greater, and it has a significant endowment of gold.

Four generalized hydrogeologic cross sections were constructed along transects in northern Nevada to test the hypotheses discussed in the Introduction (Fig. 6). The cross sections were based on an integration of existing surface and subsurface geologic and geophysical data, including geologic maps (Stewart and Carlson, 1976; John and Wrucke, 2003; Moore, 2002), published cross sections (Struhsacker, 1980; Layman, 1984; Moore and Norby, 2002; Hickey et al., 2003a, 2003b; John et al., 2003), well data (Hess et al., 2004), and models of gravity data that estimate the thickness of Cenozoic deposits (Saltus and Jachens, 1995). Stratigraphy portrayed on all of the sections was generalized into a limited number of hydrogeologic units composed of stratigraphic units inferred to have similar hydrologic properties. Geologic structure was similarly generalized such that the only structures portrayed on the sections are either those that have sufficient offset to juxtapose different hydrogeologic units or structures that were explicitly included as part of the flow scenario being tested. In certain cases, structures that were subparallel to the section were omitted or generalized as a fault more nearly perpendicular to the section trace. Although based on surface and subsurface geologic and geophysical data, the four hydrogeologic cross sections are intended to represent generalized scenarios that contain the salient geologic features to test a flow hypothesis; they are not intended to be strictly accurate representations of subsurface geology in all cases. The representation of high- and low-permeability fault zones was based on evidence from other studies and evolved as we constructed models of different scenarios. We selected the depths of the fault zones based on the need to bring hot fluids >200 °C up to the surface. Models that included fault zone depths of <5 km were “nonstarters” (not presented here), assuming a representative geothermal gradient for the Northern Nevada Rift (Blackwell, 1983). Faults zones that cut clay-rich siliciclastic units were often assigned low permeability values because this fits conceptual models of fault permeability (Bense and Person, 2006) and because of the absence of surface sinter deposits at Carlin. Although we did not represent Eocene plutons in the Carlin section, we discuss the influence of magmatic heat flow on the models. We considered the effects of Miocene dikes for the Mule Canyon section by imposing localized enhanced basal heat flow near the Muleshoe fault zone. We did this because preliminary amagmatic models of Mule Canyon were not able to match observed fluid inclusion temperature estimates.

A regional north-northwest–trending section (section A–A′, Figs. 1 and 6) portrays generalized subsurface geology along ~120 km from the general vicinity of Eureka, Nevada, at its southeastern end to Crescent Valley, near Battle Mountain and Beowawe, at its northwestern end. This section was constructed to determine whether a regional topography-driven flow system could explain the Battle Mountain heat flow high and Eureka heat flow low anomalies. The section extends to a depth of ~12 km so that deep regional groundwater flow through the Paleozoic carbonate rocks could be modeled. The section includes the following hydrogeologic units: a deep, low-permeability unit that includes Proterozoic crystalline basement rocks and Jurassic granitic plutonic rocks (unit 1, Fig. 6); low permeability Late Proterozoic–Early Cambrian siliciclastic rocks (unit 3, Fig. 6); Cambrian–Devonian dominantly carbonate rocks (unit 4, Fig. 6); low-permeability strata of the Roberts Mountains allochthon (unit 6, Fig. 6); Miocene–Pliocene sedimentary and volcanic rocks (unit 7, Fig. 6); and recent alluvium (unit 8, Fig. 6). The section has been constructed subparallel to, and just east of, the Northern Nevada Rift (Zoback et al., 1994; John et al., 2003). The section is nearly perpendicular to major range-front normal faults on the northern flank of the Roberts Mountains and Cortez Range, as well as large faults within Pine Valley and Crescent Valley (Fig. 1). The Paleozoic carbonate rocks are inferred to be reasonably continuous from the southeast end of the section near Eureka to near the center of the section in Pine Valley. The carbonate rocks are disrupted in the northern half of the section by major faults and by the Jurassic pluton that underlies the Cortez Range.

A generally east-west cross section (B–B′, Figs. 1 and 6) was constructed across the northern end of the Shoshone Range to portray generalized sub-surface geology in the vicinity of the middle Miocene Mule Canyon deposit (John et al., 2003) and the active geothermal system at Beowawe (Zoback, 1979). The purpose of this cross section was to assess the role of faults in controlling the transient hydrothermal conditions observed at Beowawe. Section B-B′ includes the following hydrogeologic units: a deep, low-permeability unit that includes Proterozoic crystalline basement rocks (unit 1, Fig. 6); low-permeability Late Proterozoic–Early Cambrian siliciclastic rocks (unit 2, Fig. 6) that include a middle 1.5-km-thick interval of higher permeability, highly fractured quartz sandstone (unit 3, Fig. 6); Cambrian–Devonian dominantly carbonate rocks (unit 4, Fig. 6), low-permeability strata of the Roberts Mountains allochthon (unit 6, Fig. 6); Miocene volcanic rocks (unit 7, Fig. 7); and alluvial basin fill (unit 8, Fig. 6). A 100-m-thick interval of potentially enhanced horizontal permeability was added at the top of the Paleozoic carbonate unit (unit 5B, Fig. 6). This interval represents a sequence boundary that commonly contains karst features and is highly permeable. Such sequence boundaries are present at various stratigraphic levels within the Paleozoic section (Cook and Corboy, 2004); the boundary was placed at the top of the hydrogeologic unit because many Carlin-type gold deposits occur at this position. We chose to put this unit at the top of unit 4, but it could have been placed within that unit. In this section, fault zones have also been distinguished as hydrologic entities based on the units they cut and the presumed degree of development of fault-related damage zones. In general, faults that cut strata of the Roberts Mountains allochthon (unit 6, Fig. 6) and Cambrian–Devonian dominantly carbonate rocks (unit 4, Fig. 6) are assigned higher fault zone permeabilities (unit 9B, Fig. 6) than fault zones developed at depth within the Late Proterozoic–Early Cambrian siliciclastic rocks (unit 10B, Fig. 6). However, in the Eocene Carlin model, faults cutting the Roberts Mountains allochthon were assigned lower permeabilities because hydrothermal fluids did not ascend into the allochthon along them.

The inferred subsurface structural geometry portrayed on section B–B′ is generalized from sections shown in John and Wrucke (2003) and in Tilden et al. (2005). Major structural elements include west-dipping range-bounding faults on the west side of the Shoshone Range, steep faults in the vicinity of the Mule Canyon deposit that proxy for the east-dipping Muleshoe fault (John et al., 2003) and for highly anisotropic permeability within dike swarms of the Northern Nevada Rift (John and Wrucke, 2003). A single steeply west dipping fault is portrayed in the vicinity of Beowawe; this fault proxies for the more gently dipping Malpais fault and generalizes the east-northeast–trending and north-northwest–trending faults that localize geothermal fluids (Zoback, 1979; Struhsacker, 1980). All of these faults are portrayed as cutting both low-permeability strata of the Roberts Mountains allochthon (unit 6, Fig. 6) and underlying Cambrian–Devonian dominantly carbonate rocks (unit 4, Fig. 6). The top of the Paleozoic carbonate rock section is present at maximum depths of ~6 km below land surface.

Section C–C′ (Figs. 1 and 6) is located along the same line of section as cross-section B–B′. This section portrays the inferred geology of the northern end of the Shoshone Range, including the Mule Canyon deposit (John et al., 2003) as it might have appeared in middle Miocene time, during Mule Canyon mineralization. The purpose of this section was to compare fluid circulation from the present-day Beowawe scenario with fluid circulation at the time of low-sulfidation Au-Ag ore deposition at Mule Canyon. Stratigraphic and structural elements in this section are generalized after those shown in John et al. (2003, their Fig. 18). Hydrogeologic units are the same as those described for section B–B′. Restoration of section B–B′ to middle Miocene time involved removing range-front fault offset on the west side of the Shoshone Range; this area is shown on the west side of C–C′ as an uplifted block of Paleozoic rocks, although there may have been a cover of volcanic rocks at that time. The restoration also removed offset on younger, east-northeast–trending faults east of the Muleshoe fault, reducing the amount of fault offset present near Beowawe. While displacement along the Malpais Rim is late Miocene in age (John et al., 2003), our hydrostratigraphic reconstruction includes a range- bounding fault as a permeable conduit for fluids. Also shown on section C–C′ is a shallow pluvial lake, as suggested by John et al. (2003). The middle Miocene restoration has the effect of reducing fault offset and making the Paleozoic carbonate rocks more continuous at depth. An important difference between the Mule Canyon and Beowawe sections is that the distance between framing faults is much farther apart along C–C′relative to B–B′. This is important in controlling the fluid-rock isotope composition of the discharging fluids.

A fourth cross section (D–D′ in Figs. 1 and 6) was constructed along a line trending northeast-southwest that is intended to represent the hydrogeologic setting of the Carlin trend during Eocene time. The purpose of this cross section was to see if the fault-controlled fluid flow system similar to what is found today at Beowawe and during Miocene time at Mule Canyon could have formed the world-class Carlin ore district. Our section does not include plutonic intrusions that are inferred to underlie the Carlin trend at a depth of between 5 and 10 km based on aeromagnetic features (Ressel and Henry, 2006). We used 10 hydrostratigraphic units to represent the siliciclastic, carbonate, and fault units. The lowermost unit (unit 1) is made up of Proterozoic metamorphic rock, which is overlain by Proterozoic and Paleozoic lower plate siliciclastic rocks. The siliciclastic rocks include low-permeability interbedded shale (unit 2) and sandstone that surround an interval in the middle (unit 3), a 1.5-km-thick fractured quartzite unit with potential fracture permeability. This fractured interval is intended to provide a deep potential fluid pathway for ore-bearing flu-ids. Unit 4 is made up of lower Paleozoic carbonate rocks. As with the Beowawe model, this unit is assumed to have a moderate to low permeability, we have also included a thin permeable karst interval (unit 5) at the top of unit 4. This interval is quite thin, 40 m or less. Its position within the carbonate unit is arbitrary. It could also represent a subhorizontal fault zone. A key difference between the Carlin and Beowawe sections is that the carbonate unit is shallow, and is therefore more likely to receive warm upwelling groundwater. On the Beowawe–Mule Canyon sections, the carbonate unit is deeper and (owing to its position) is not part of the discharge area. Three steeply dipping, permeable normal faults cut the entire succession (unit 9). Faults are shown as planar features that dip ~80° east and west. Where the faults cut a low-permeability rock unit, the permeability of that fault segment is believed to be low (unit 10). The central fault is intended to represent the Post-Gen fault zone.

Mathematical modeling represents an important method for providing insights into crustal-scale fluid-flow patterns during metallogenesis (Garven et al., 1993; Person et al., 1996). We solved equations representing groundwater flow, heat, and silica transport as well as groundwater residence times (Goode, 1996) and fluid-rock isotope exchange using a parallel version of the model described by Person et al. (2007). For silica transport, we incorporated the kinetics of quartz precipitation following the approach of Rimstidt and Barnes (1980). We also computed apparent apatite fission track ages along a series of points near the Carlin trend to determine if this thermochronometer could be used to constrain paleo-hydrologic model results for this system. The transport equations were solved using the finite element method. Finite element models were constructed for each of the cross sections presented in Figure 6. The geologic cross sections were discretized using between 1116 and 3522 nodes. The length of the triangular elements in the x dimension varied from ~1000 m in the low-permeability crustal rocks to 40 m in the fault zone. Time step size varied between 0.01 and 10 yr for the different cross sections. For the Eureka–Crescent Valley (A–A′) cross section, we experimented with the permeability and anisotropy of each hydrostratigraphic unit in an attempt to match modern heat flow constraints described in the discussion of the hydrogeology of the Great Basin. For the Beowawe, Mule Canyon, and Carlin cross sections (B–B′, C–C′, and D–D′), in addition to varying sediment permeability and permeability anisotropy, we also varied the permeability and anisotropy of fault units (9–10). We found that we could best match modern heat flow data by assigning different permeabilities to segments of the subvertical fault zones. This was done to qualitatively relate fault permeability to the type of rocks cut by the faults (Bense and Person, 2006). We found that to match modern thermal and isotope data, faults cutting fine-grained siliciclastic units needed to be assigned relatively lower permeability values (three orders of magnitude) than faults cutting carbonate or fractured siliciclastic units. We then used these values to inform our selection of fault and sediment permeability values in our paleohydrologic reconstructions of the Miocene and Eocene gold mineralization models. We found that using similar hydrogeologic framework models for both modern and fossil flow systems increased our confidence in our findings. We treated faults as conduit-barrier systems (Bense and Person, 2006). We assigned maximum permeability in the direction of the fault plane (k_max; see equation A3 in  Appendix) and low permeability (k_min) perpendicular to the direction of the fault plane.

To satisfy the grid courant number, each simulation had to be run for tens of thousands of time steps. Given the large number of computational time steps for solute transport, and multiple fluid and/or rock tracers (SiO₂, groundwater residence times, heat, and δ¹⁸O), our parallel computing technique reduced the overall computational time by a factor of four (Person et al., 2006) using an eight processor Pentium-4 Linux server. The governing flow and transport equations used in this study are presented in the  Appendix. Below we summarize the rock properties and isotopic input data, and silica dissolution and/or precipitation kinetics used in our model.

The reactive surface areas used were based on fracture aperture density rather than mineral grain size (Brantley and Mellott, 2000). We set an aperture density of 1 m for fault zones and 0.1 m for other hydrostratigraphic units, respectively which also defines the reactive surface area for isotopic and silica fluid-rock reactions. While these values are arbitrary, setting a surface area based on grain size would yield fluid-rock isotope exchange rates far higher than observed.

The thermal and hydrologic properties used in this model for the crystalline and sedimentary rocks are listed in Tables 5 and 6, respectively. The permeability values were largely based on model calibration, but are consistent with values reported in the literature for fractures and sedimentary rocks. The thermal properties used are representative of sedimentary, volcanic, and metamorphic rocks (Person et al., 1996). Excluding coals and evaporates, the thermal conductivity of sedimentary rocks varies between ~1.0 and 3 W m⁻¹ °C⁻¹ (Nunn and Lin, 2002). For most sedimentary rocks, and thermal conductivity variations are largely due to porosity variations rather than to absolute variations in thermal conductivity of the mineral phases. The heat capacity of solid and fluid phases as well as fluid thermal conductivity varies modestly with temperature. We chose not to represent temperature effects on heat capacity variations because they have a second-order effect on computed subsurface temperatures (Table 5). Dispersivity is known to vary widely with scale (Gelhar, 1993). Because our focus was on fault-controlled flow and transport, we selected low values of longitudinal and transverse dispersivity (Table 5). Fault zone permeability can be influenced by a number of factors, including the development of a low-permeability core zone composed of fault gouge (Morrow et al., 1984), a higher permeability outer damage zone (Caine et al., 1996), syndeformational mineralization, and the regional stress fields (Hobbs et al., 2004). Within sedimentary rocks, smearing of clay and sand facies into fault zones can result in faults developing a conduit and/or barrier behavior with the highest permeability parallel to the fault plane and low permeability normal to the fault plane (Bense and Person, 2006). We used published permeability data sets for crustal rocks as a guide in selecting these values (Clauser, 1992; Brace, 1980, 1984). We varied the permeability of crustal rocks between 10⁻¹⁷ to 10⁻¹⁹ m² outside the fault zone. Fault zone permeability (k_z) was varied between 10⁻¹⁴ and 10⁻¹¹ m².

Fault permeability for regions of the cross-sections B–B′, C–C′, and D–D′ that cut carbonate rocks were assigned a vertical permeability of 10⁻¹¹ m². Regions of fault zones cut by fine-grained siliciclastic sediments were typically assigned a permeability of 10⁻¹³ m². For the Carlin cross section, subvertical faults were assigned a vertical permeability of 10⁻¹³ m². Where subvertical faults cut the Roberts Mountains allochthon, they were assigned a permeability of 10⁻¹⁵ m². The permeability of the fault in the horizontal direction was assigned a value 10 times less than the vertical direction.

We used a representative value of permeability (10⁻¹⁷ m²) to model flow within crustal igneous and metamorphic rocks. Recent permeability measurements obtained from deep boreholes (Stober and Bucher, 2005) as wells as regional syntheses (Manning and Ingebritsen, 1999) suggest that crystalline rock permeability can be as high as 10⁻¹⁵ m² or higher in the upper few kilometers of the Earth's crust. Using a higher permeability value for crustal rocks in our model would have enhanced the potential for deep convective heat transfer below the sedimentary pile. However, it is important to remember that the presence of overlying low-permeability siliciclastic confining units would prevent crustal fluids from rising to the surface without the presence of permeable fault conduits. In addition, the fault permeability values we used were orders of magnitude (10⁻¹¹ to 10⁻¹³ m²) higher than the highest estimates for bulk crustal permeability reported by these studies.

All subvertical faults are relatively thin (~40 m) to match observed thermal, isotopic, and geochemical features; subvertical faults for the Beowawe models are connected hydrologically to a thin aquifer unit at the top of the Paleozoic carbonates (locally karst), and, in the Carlin trend, by a deep, thick, more porous aquifer in the middle of the siliciclastic sequence (10⁻¹³ m²). We assigned a fault zone permeability anisotropy of 100 (k_z/k_x) to match observed temperatures and the narrow width of geothermal anomalies. Simulations that incorporated an isotropic fault zone were unable to produce the narrow high heat flow zones observed or isotopically light oxygen isotopic composition. The karst unit was assigned an anisotropy of 10 (k_z/k_x) so as to facilitate lateral flow. Practical considerations regarding simulation time restricted our ability to represent the fault zones using grid refinement on the submeter scale, which would have been desirable.

For groundwater flow we assumed that relief on the water table is a subdued replica of the land surface topography (Toth, 1963). We specified the head at the top surface for all three cross sections. No-flow boundaries were assumed along the base and sides of the solution domain. A no-flow boundary was applied to the bottom and side boundaries. For stream functions, we specified a constant value (zero) stream function boundary along the sides and base of the domain. We specified a flux (dΨ/dz) across the top boundary proportional to the lateral water-table gradient as described by Senger and Fogg (1990a, 1990b). For heat flow, a specified temperature of 10 °C was imposed along the water table except where fault zones occurred and groundwater flow was upward. Here a convective boundary condition (i.e., ∂Τ/∂z = 0) was imposed. This was done for the Malpais fault along cross-section B–B′, the Muleshoe fault along C–C′, and the Post-Gen fault along D–D′. For B–B′, C–C′, and D–D′, a basal heat flux of 104 mW/m² was imposed on the base of the domain consistent with the Battle Mountain heat flow high. Along the Muleshoe fault along C–C′, heat flow decayed exponentially from 300 to 104 mW/m² over a period of 20 k.y. to represent the cooling of Miocene mafic dikes. For A–A′, a basal heat flux of 60 mW/m² was imposed. We allowed the initial fluid isotopic composition to be in equilibrium with the rocks assuming a conductive temperature field. We allowed fluid to infiltrate with an isotopic composition of –14‰ for the present-day Beowawe and Miocene Mule Canyon model (Friedman et al., 2002) and –18‰ for the Eocene Carlin cross section (D–D′) when temperatures in this region are thought to be cooler (Norris et al., 1996). Estimates of Holocene–Miocene meteoric precipitation are also consistent with data from the ore deposits (Table 1). An advective boundary condition (dR¹⁸O/dz = 0) was imposed at the water table along faults that served as discharge areas. For silica transport, we set the concentration at the water table to zero and imposed no flux conditions along all other sides. While chemical weather produces some dissolved silica in shallow groundwaters (Feth et al., 1964), these values are low when compared to geothermal reservoir conditions (Rimstidt and Barnes, 1980; Rimstidt and Cole, 1983). For the fault zones that were discharge areas, an advective boundary condition was used for silica transport as described above.

We first present hydrothermal model results to test the hypothesis that regional groundwater flow within the Paleozoic carbonate aquifer system is responsible for the observed Battle Mountain heat flow high (>104 mW/m²; BMH in Fig. 2) and the Eureka heat flow low (<60 mW/m²; EL in Fig. 2). It is important to evaluate this hypothesis because regional groundwater flow systems capable of producing thermal anomalies may generate world-class ore deposits (Garven et al., 1993). The hydrologic cross-section A–A′ in Figure 7 is more or less perpendicular to the regional water-table gradient across the Upper Humboldt River regional flow system (see head contours, Fig. 2). The cross section extends a few kilometers into the recharge area of the Railroad Valley regional flow system (Fig. 2).

We intentionally did not incorporate any permeable faults or karst layers along A–A′ to focus on the effects of the regional, topography-driven flow. The most permeable unit represented in our model was the Paleozoic carbonate aquifer (unit 4A; 10⁻¹⁴ m²). While aquifer tests report much higher permeabilities for some members of the Paleozoic carbonate aquifer, our model lumps these higher permeability units together with lower permeability members. We imposed an anisotropy of 1000 (kx/kz) for the carbonate unit that combines the effects of high-permeability karst layers embedded in lower permeability carbonates. Had we increased the vertical carbonate permeability, the overall match to observed heat flow conditions would not have been good (i.e., there would have been an overall bias in the average surface conductive heat flow).

For the A–A′ model (Fig. 7), we only present results for groundwater flow and heat transport. A series of local flow cells, as indicated by stream functions in Figure 7C, formed in the southern portion of A–A′ where the Paleozoic carbonate rocks crop out (Fig. 7). Water-table topographic relief in this area resulted in the formation of local recharge and/or discharge areas. No free convection cells formed across this cross section. This produced a series of shallow heat flow anomalies (dashed line in Fig. 7A) that ranged between <40 and 100 mW/m² between recharge and discharge areas, respectively. While the magnitude of these anomalies is similar to observed levels in the Eureka heat flow low and Battle Mountain heat flow high (dashed red lines in Fig. 2), the patterns are largely inconsistent. We conclude that the regional Battle Mountain heat flow high anomaly cannot be explained by convective heating. We note, however, that the areas of low computed heat flow correspond to regions where Paleozoic carbonates crop out in the south and are largely consistent with the position of the Eureka heat flow low in Figure 2.

Prior studies of the Beowawe (Faulder et al., 1997) and Dixie Valley (McKenna and Blackwell, 2004) geothermal systems concluded that these fields are the result of focused, transient, fault-controlled groundwater circulation. The formation of sinter deposits (source reservoir temperatures >180 °C) indicates that fluid circulation must be very deep (Rimstidt and Cole, 1983). Groundwaters within the discharge area show very little isotopic enrichment, suggesting low surface areas, low Damkholer numbers, and high velocities (Person et al., 2007). Our models test the hypothesis that subvertical faults, which result from extensional or transtensional stresses, can serve as conduits for groundwater, heat, and dissolved silica (e.g., see López and Smith, 1995; Sibson et al., 1988; Sibson and Scott, 1998).

We constructed hydrothermal models along B–B′ to (1) assess the effective permeability and permeability anisotropy of the Malpais fault zone; (2) determine the depth of fluid circulation required to form shallow heat flow anomalies in excess of 2000 mW/m²; (3) constrain the duration of the flow system; and (4) determine whether spatial variations in basal heat flow are needed to account for such a focused shallow heat flow anomaly. Because faults near Beowawe were generalized as a single steeply dipping structure, the thin (40 m) high-permeability karst unit was critical for connecting adjacent subvertical fault zones to create a single pass convection cell. Modeling a gently dipping Malpais fault would have reduced the role of the thin karst zone by allowing the Malpais fault to intersect the Muleshoe fault at depth, but this geometry would not dramatically change the shape or size of the modeled convection cell. We present one of a suite of simulations considered by Banerjee et al. (2007). The simulation was run for 20 k.y. using a fluid flow time-step size of 0.1 yr. Because we used a Lagrangian particle tracking scheme (modified method of characteristics) to solve for both enthalpy and solute transport, 10–200 particle move time steps were required for every fluid flow time step. We used subvertical fault zone permeability similar to what was measured by hydraulic testing (10⁻¹¹ m²).

Computed stream functions after 6 k.y. (Fig. 8A) indicate that almost all of the flow up the Malpais fault originates along a subvertical fault to the west along the Argenta Rim (i.e., down the Muleshoe fault). Stream functions vary from ~100 to –100 m²/yr. Negative stream functions denote flow from left to right while positive stream function values denote flow from right to left. It is important to point out that the lateral circulation distances at Beowawe are relatively short (~5 km), providing relatively little opportunity for fluid-rock isotope exchange. Owing to the high vertical permeability within the fault plane, groundwater flow rates ranged between 30 and 140 m/yr during the simulation within the Malpais fault (Fig. 9B). The depth of circulation is constrained by the position of a thin, permeable karst horizon (unit 5B, Fig. 6). Groundwater flow rates were highest at the start of the simulation (Fig. 9B). Decay in convective heat flow occurs after the initial fluids within the fault zone are replaced by cooler meteoric fluids from the Argenta Rim that descend along the Muleshoe fault (Fig. 9A). By 5 k.y., fluids along the Malpais fault are nearly isothermal at depth (Fig. 10A). Comparison of computed and observed vertical temperature profiles for this simulation at different distances to the east of the Malpais fault after 300 yr have similar magnitudes and gradients, although there is not a perfect match (Fig. 4B). This is likely due to three-dimensional nature of flow along the Malpais fault zone. The silica and isotopic composition of fluids within most hydrostratigraphic units (Figs. 8B, 8C) reflects the initial equilibrium conditions with the host rock at a given temperature. Retrograde solubility of silica above temperatures of ~300 °C accounts for the decrease in dissolved silica below a depth of ~9 km. Fluids within the Malpais fault zone show strong advective affects (Fig. 8C) due to vigorous advective transport. Fault zone pore fluids are not in equilibrium with the surrounding wall rocks. Fault zone fluids become depleted in δ¹⁸O, but only by ~8‰ (−7‰) near the surface after 20 k.y. (Fig. 9E). The initially enriched (+8‰) isotopic composition of the discharging fluids at early time reflects the initial isotopic compositions of higher temperature reservoir fluids (which we assume to be in thermal and isotopic equilibrium with the surrounding mineral phases at the beginning of the simulation). A time of ~5 k.y. is required to replace these reservoir fluids in the fault and/or karst zone with relatively less enriched meteoric fluids. This is more enriched than observed and is probably due to grid dispersion along fault elements. Grid dispersion occurs when the grid size exceeds the width of the solute front, which is in turn controlled by the dispersivity. It is generally desirable to have grid refinement comparable to the magnitude of the longitudinal dispersivity (in our case, ~1 m) along the fault zones where advective transport dominates. Because of computational limitations, our fault element size was ~40 m. High silica concentrations and isotopically enriched fluids along western and easternmost faults are induced by density gradients and suggest that mineralization could occur far outside the central rift zone.

Dissolved silica concentrations within the Malpais fault are initially high (~800 ppm or 0.013 mol/m³) at the start of the simulation, reflecting the initial assumption of geochemical equilibrium with the rocks (Fig. 9C). These are flushed out by cooler fluids with lower dissolved silica concentrations. After 5 k.y. concentrations drop to ~300 mg/L, reflecting the equilibrium solubility of quartz at temperatures above 200 °C (Figs. 9C and 10B). After 20 k.y., silica concentration is ~205 mg/L. Kinetic-based quartz precipitation could not keep pace with advective transport (i.e., low Damkholer numbers). Silica concentrations increased across the high-permeability carbonate aquifer because undersaturated, relatively cool fluids moving down from the Argenta Rim were heated as they moved laterally. Peak silica concentrations within the Malpais fault occurred between 100 and 300 yr after the onset of convection (Figs. 9D and 10B). Given enough time, silica concentrations within the fault zone will decline due to the high flow rates along the fracture planes. Computed fluid temperatures along the Malpais fault zone near the water table decreased from ~170 to 110 °C between 300 and 5000 yr (Fig. 10A). Temperatures become increasingly isothermal as the simulation progresses.

Assuming a fault width of 40 m and a fault length of 3.75 km (Rimstidt and Cole, 1983), after 20 k.y. ~1.3 × 10¹⁰ kg of silica have been brought to the surface (Fig. 9D) and would potentially be available for the formation of sinter deposits. This is about one order of magnitude lower than what is observed at Beowawe (Rimstidt and Cole, 1983). This suggests that the flow system that formed the sinter deposit occurred over a longer period of time (~200 k.y.). However, it is unlikely that flow was continuous because simulated temperatures and silica concentrations would decline as meteoric fluids cooled the system. If our model is accurate, multiple episodes of fluid flow were required to form the Beowawe sinter deposits. The Neotectonic analysis of Wesnousky et al. (2005) shows an earthquake recurrence interval along the Malpais fault of 10 k.y., consistent with our hypotheses assuming that tectonic activity is responsible for generating fault permeability and fluid flow (Micklethwaite and Cox, 2004, 2006). These episodic flows would likely generate more shallow silica mineral precipitation than continuous flow because intervening periods of quiescence would allow pore fluids to equilibrate at higher temperatures (due to convective cooling) with silica-bearing minerals at higher temperatures. Thus, shallow discharge would on average have a higher dissolved silica concentration during the first few thousand years, as shown in Figure 9E. It is also possible that multiple fault zones (that strike at high angles to the Malpais fault) are involved in sinter formation at Beowawe.

Convective temperature anomalies form within a few kilometers of the fault zone after ~6 k.y. (Fig. 11C). A similar area of convective cooling is located next to the feeder fault near the Argenta Rim. This is largely consistent with field observations that show that no area of low heat flow exists within Whirlwind Valley within 10 km of the Malpais fault zone. It is important to note the significant degree of convective cooling along this single pass flow cell depicted in Figure 11. After 20 k.y. (not shown), at least two-thirds of the flow cell has become cooler than initial conductive conditions. This suppresses the amount of fluid rock-isotope exchange and dissolved silica concentrations.

The goal of this simulation was to determine what combination of Miocene hydrologic and/or thermal boundary conditions and fault and/or sediment permeability would be consistent with the extremely high temperatures (~200 °C at 150 m depth) and enriched oxygen isotopic composition (10‰ shift) along the Muleshoe fault at a relatively shallow depths. We were also interested in determining whether the same hydrologic framework model could account for both modern and fossil heat flow conditions. We used the same fault permeability values for the subvertical faults (10⁻¹¹ m²). For this simulation (C–C′), we had to enforce a lake boundary condition (i.e., constant head of 3,500 m across the valley) so that upwelling occurred along the Muleshoe fault (see Fig. 6). Recharge was inferred to have occurred along a range-front fault ~15 km east of Mule Canyon during Miocene time. We also found that to approach the maximum fluid inclusion homogenization temperatures measured at Mule Canyon, we had to impose a spatially varying basal heat flux boundary condition (consistent with dike emplacement) in the vicinity of the Muleshoe fault. Heat flow was locally three times higher (300 mW/m²) than background values, but allowed to decay over a period of 20 k.y. This produced temperatures along the Muleshoe fault of ~150 °C at 280 m depth (shallowest node below the top surface along Muleshoe fault). Temperatures reached a maximum within ~5 k.y. and then stabilized (Fig. 12A). Even higher imposed basal heat flux (perhaps 600 mW/m²) would be needed to exactly match the shallow fluid inclusion temperatures. Maximum computed basal temperatures along the Muleshoe fault were ~750 °C. It is important to note that the high heat flow was not needed to change the sense of circulation. However, the high heat flow was needed to approach the observed fluid inclusion temperatures. Evidence that the Muleshoe fault became sealed during mineralization is not accounted for in this model (John et al., 2003). This led to higher temperatures at shallow depths as fluid pressures approached lithostatic conditions. Significant cooling occurred in the vicinity of the eastern range-front fault (Fig. 12B). Computed fluid oxygen isotopic compositions were more shifted (~20‰; Fig. 12D) than the fluids that formed the deposits (~10‰; Table 1). The greater shift in fluid isotopic composition between this model and the modern Beowawe system is likely due to the longer flow path (17 km compared to 5 km) of fluids across this Miocene flow system (Fig. 12C). Computed fluids near the lake sediment-water interface had high silica concentrations (not shown; ~500 ppm), but high temperatures probably prevented much silica precipitation.

Here we consider hydrothermal flow in the Eocene Carlin trend along section D-D′, and evaluate whether the hydrologic framework model presented above for the Beowawe and Mule Canyon systems can account for the isotopic data, fluid inclusion temperatures, and apatite fission track annealing found at Carlin. Banerjee et al. (2007) found that while the broad features of the Beowawe model were consistent with Carlin thermal and isotopic constraints, the vertical fault permeability had to be decreased by at least two orders of magnitude (10⁻¹³ m²) to attain the required temperatures. Because the Post-Gen fault intersects low-permeability sediments of the Roberts Mountains allochthon near the surface, we lowered the fault permeability as it cut this siliciclastic unit by two orders of magnitude (D–D′; Fig. 6). This, in effect, resulted in the formation of a series of closed-loop natural convection cells. In this model, fluid circulation extends far below the carbonate layer into permeable Neoproterozoic siliciclastic rocks. These sedimentary rocks are a possible source of gold for the ore-forming fluids (cf. Ilchik and Barton, 1997).

Computed stream functions (Fig. 13) indicate that flow patterns evolved on time scales of 10⁴ yr due to changes in fluid density associated with convective heat transfer effects. Local convection cells within permeable siliciclastic layers formed during the first 30 k.y. These local cells had length scales of ~1–2 km within the siliciclastic rocks and occurred both at relatively shallow and intermediate crustal depths (2–8 km). A single more laterally and vertically extensive (12 × 6 km) convection cell also formed during the first 30 k.y. on the eastern side of D-D′. Convection in this larger cell occurs within subvertical faults, the thin (40 m thick), relatively shallow karst unit, and the deeper, relatively porous, high- permeability Neoproterozoic siliciclastic unit. The depth of fluid circulation was 6–8 km in this cell. By 40 k.y. a second large cell formed in the western portion of D–D′ at the expense of the smaller convection cells. These large cells gained strength in response to temperature redistribution (cf. Figs. 13A and 13C). After 40 k.y., the two large convection cells dominated. The maximum vertical velocity in the upflow zone of the subvertical fault was 7.8 m/yr. It is important to point out that the lateral scales of these flow cells are about twice the lateral length scale of the Beowawe single pass convection cell. It is also important to note that the Carlin trend convection cells are more or less closed loops rather than single pass cells. Only modest amounts of meteoric water were added to the convection loop along the western subvertical fault.

The thermal evolution of D–D′ showed significant convective anomalies developing along the subvertical faults (Fig. 14). By 30 k.y., temperatures reached 250 °C at a depth of 500 m within the Post-Gen fault. Shallow conductive surface heat flow was ~600 mW/m². The footprint of this positive heat flow anomaly extended laterally for ~5 km during the first 30 k.y. At later times, a broad, negative heat flow anomaly (i.e., <60 mW/m²) formed as a result of downward flow along the subvertical faults east and west of the Post-Gen fault. By 200 k.y., lateral conductive cooling associated with descending cool meteoric groundwater flow along the eastern and western subvertical faults extended the region of low heat flow to more than two-thirds of the cells lateral footprint (one third of the cell has a positive thermal anomaly). This acted to slow the rate of fluid-rock isotope exchange.

The formation of large-scale convection cells along D–D′ resulted in significant fluid-rock isotope exchange in the fault zones when temperatures exceeded 250 °C (Figs. 15A, 15B). Isotopic fluid composition within the convection loop varied between –8.6‰ and +4‰. This is due to both mixing of deep and shallow fluids as well as fluid-rock isotopic exchange. Within the upflow zone, the oxygen isotopic composition was substantially enriched relative to the Beowawe upflow zone fault system. This is because the permeability and vertical velocity within this fault zone are much lower. Because the convection cells include a region of relatively high porosity siliciclastic sedimentary rocks at depth, significant fluid mixing occurs between hot isotopically evolved fluid and cool unevolved meteoric water. Isotope exchange with the Paleozoic carbonate rocks within the discharge area occurred principally within the high-permeability carbonate aquifer as well as adjacent to the Post-Gen fault. The isotopic composition of the carbonate rocks declined from 22‰ to 8‰ over a period of 200 k.y. The magnitude and pattern of fluid-rock isotope exchange near the intersection of the Post-Gen fault with the high-permeability aquifer (Fig. 15B) mimics the range of carbonate δ¹⁸O values (Hofstra and Cline, 2000) and alteration patterns observed in the deposits (Fig. 5A). Computed groundwater residence times are presented in Figure 15C. What is remarkable about the computed age distribution is that relatively young meteoric fluids do not dominate within the regional convection cells. This suggests that only a minor amount of unexchanged meteoric water was introduced to the convection cells. The meteoric water is added along the eastern and western bounding faults of the section to the east and west of the Post-Gen fault, as evidenced by the relatively young water and isotopically depleted fluids. In other words, this is not a single pass system. When this system was active (39 Ma), the age of the groundwater within the convection cells was ~190 k.y. older, approaching the maximum possible age for groundwater in this model of 200 k.y. Note that we used an exponential scale in contouring groundwater age. This result stems from the low permeability assigned to the Post-Gen fault zone that cuts the fine-grained Roberts Mountains allochthon above the carbonate unit.

We monitored the time-temperature history of four nodes just above the Paleozoic karst unit along D–D′ (points a–d shown on cross-section D–D′; Fig. 6). These points spanned the groundwater recharge to discharge area on the western convection cell. Changes in temperature during the 200 k.y. hydrothermal event were superimposed onto the much longer (114 m.y.) conductive cooling and/or uplift related thermal history presented by Cline et al. (2005) for the Carlin trend (Fig. 16B). We used equations 10–13 to compute the apparent fission track age and fission track distribution for points a–d (Fig. 16A). Computed fission track ages varied between 91.6 and 35.5 Ma between groundwater recharge and discharge areas, respectively (Fig. 16A). We interpret modeled fission track ages of 35.3 Ma as indicating complete resetting of fission tracks older than the age of mineralization in areas where rock temperatures exceeded 140 °C. It is important to note that annealing continued after the 39 Ma age of mineralization as rock temperatures conductively cooled as they moved toward the land surface between 39 and 0 Ma. At a distance of 4 km away from the Carlin fault, a bimodal fission track distribution was generated due to convective thermal effects and partial annealing of fission tracks older than the age of mineralization (points b and c). The lack of spread of track lengths for the recharge area (point a) is likely due to the relatively isothermal conditions between 60 and 20 Ma associated with tectonic stability. The simulated apparent fission track age of 35.3 Ma is younger than the mineralization event due to the continued conductive cooling and annealing modeled between 39 and 0 Ma.

This paper has summarized our efforts to understand the plumbing of active geothermal systems: the Miocene volcanic-hosted Mule Canyon Au-Ag deposits, and the Eocene carbonate-hosted, Carlin-type deposits. We used a “present is the key to the past” approach in that we first developed a conceptual hydrogeologic framework model of modern Great Basin hydrothermal systems (Eureka–Crescent Valley, A–A′ and Beowawe, B–B′). We tested this model using fluid enthalpy, isotopic, and chemistry data. We then applied this hydrogeologic framework model to reconstruct past (Miocene and Eocene) hydrologic conditions. Using this approach, we realized that to match both modern geothermal conditions and paleoflow indicators such as fluid inclusion temperatures and fluid-rock isotope exchange, faults must be highly anisotropic with higher permeability within the fault plane (and lower permeability perpendicular to the fault plane). Faults that cut fine-grained siliciclastic units needed to be assigned lower permeability values than those cutting carbonates or coarse-grained sand units. These findings are largely consistent with the fault permeability model proposed by Bense and Person (2006). We also realized that modern and probably fossil hydrothermal systems must be transient in nature. It is likely that Cenozoic–Quaternary tectonics are ultimately responsible for creating the high-permeability fault systems associated with the modern geothermal features and gold mineralization across the Basin and Range (Micklethwaite and Cox, 2004, 2006). This is consistent with the neotectonic study of Wesnousky et al. (2005), who noted an earthquake recurrence of 10 k.y.

Below we summarize in more detail our findings from each cross-sectional model run. We used a combination of measured rock permeability and model calibration (i.e., matching both inferred and measured temperatures and silica concentrations) to arrive at the set of permeability values. Likewise, rock surface area was estimated based on observed fluid and rock isotopic alteration measurements. Thus, our models must be viewed as nonunique. The ranges of permeability values we used in our model are broadly consistent with those measured from aquifer tests within the Great Basin as well as for a given lithology type (Freeze and Cherry, 1979). That said, there is a great range of uncertainty in rock permeability for a given lithology (3–5 orders of magnitude). It is likely that the permeability of the coarse-grained siliciclastic and carbonate units could have been higher by one to two orders of magnitude and still produced similar hydrothermal results. Likewise, the shale units could have had a permeability lower than what was used in our model and this would not have changed our results much. However, lowering the permeability of the sand and/or carbonate units much below the vertical permeability of the fault zones vertical permeability would have produced a bottleneck effect and lowered the flow rates within the fault zone for these cross-sectional models. The fault permeabilities used in our models (10⁻¹¹ and 10⁻¹³ m²) were within the bounds of measured fault permeabilities from aquifer tests across the Basin and Range at Yucca Mountain, Dixie Valley, and Beowawe (10⁻¹⁰ to 10⁻¹⁵ m²). It is likely that over time scales of 10³ yr, fault permeability changes due to mineralization. This was not represented in our model. While representing the fault zones using higher permeabilities would have helped match the total volume of silica contained within the sinter deposits, it would have raised the temperatures above what is observed. Thus, using multiple geochemical, thermal, and isotopic tracers provided us with more confidence in our choice of model parameters and fault plane geometry.

Results from the mathematical model along A–A' lead us to conclude that regional (~100 km) groundwater flow systems cannot be called upon to explain the present-day Battle Mountain heat flow high. We conclude that this broad anomaly must be due to crustal thinning, as concluded by Blackwell (1983) and Heimgartner et al. (2005). There is no possibility that regional groundwater flow could have produced the heat flow anomalies in Whirlwind Valley at Beowawe that exceed 2000 mW/m². On the other hand, our model results could help to explain the occurrence of the Eureka heat flow low in Figure 2 by convective cooling.

Our model results suggest that the Beowawe system represents a permeable, single-pass, relatively short (~5 km) flow system with ground-water discharge at the surface. Discharge at the Beowawe geyser field results from flow along thin, fracture-dominated flow paths. Fracture controlled flow results in higher fluid velocities, cooler temperatures, and less fluid-rock isotope exchange and hydrodynamic mixing of older, more exchanged fluids (relative to porous media flow). The longer fluid flow distances (~17 km) represented in our Miocene Mule Canyon model (C-C′) resulted in relatively more mixing and fluid-rock isotope exchange along the flow path. The reversal in flow directions during Miocene time was likely due to the presence of a shallow rift lake across Whirlwind Valley and higher heat flow associated with the emplacement of volcanic dikes near Mule Canyon.

There are important differences in the style of convective flow in the active Beowawe and Eocene Carlin systems. In the Carlin flow system, the lower permeability of the Post-Gen fault zone in the Roberts Mountains allochthon and the high permeability in underlying carbonate and siliciclastic units (D–D′ in Fig. 6) created a closed-loop convection cell. Most of the water within the flow loop recirculates and has much older calculated fluid ages than in the case of the Beowawe system. Lateral flow within porous, thick clastic rocks at the base of the system produces an isotopically exchanged fluid (ore fluid) that ascends and mixes with unexchanged meteoric water at shallow levels, as inferred from isotopic studies of the deposits.

Modern temperature anomalies found at Beowawe and Dixie Valley as well as apatite fission track data along the Carlin trend suggest that the hydrothermal events related to gold metallogenesis must have been highly transient, natural convection systems. We believe that the Beowawe system is 5 ka or younger. If the Beowawe system was longer lived, large areas of Whirlwind Valley should have anomalously low heat flow. This study identifies possible recharge areas within Dixie Valley, strongly suggesting a three-dimensional, fault-controlled, natural convection system. In the absence of active tectonics, these flow systems must shut down quickly. Mature (>30 k.y.) flow systems should have broad areas of low heat flow, which are absent everywhere within the Basin and Range. The Carlin trend system may have been longer lived due to carbonate dissolution, which created permeability. This would have occurred wherever carbonates were heated in the shallow upflowing and outflowing portions of the convection cells. The Carlin trend system may have also been sustained by added heat from magmatic activity.

The computed fission track ages (91.6–35.3 Ma) are in relatively good agreement with the distribution of apatite fission track ages collected across the Carlin trend shown in Figure 15. It is interesting to note that observed fission track ages of younger than 39 Ma occur along the Carlin trend (Fig. 15C). They suggest that younger hydrothermal events have occurred along the Carlin trend, such as the Miocene Hollister low sulfidation Au-Ag deposit (Wallace, 2003), the rhyolite flows just west of Goldstrike (15.2 Ma), and the presence of small sinter deposits in the Carlin Basin (ca. 16 Ma, Wallace et al., 2008).

With the exception of the Mule Canyon model (C–C′), our cross-sectional geologic reconstructions shown in Figure 6 did not consider the role of magmatic bodies on fluid circulation. As noted earlier, in the discussion of Miocene epithermal systems and Eocene Carlin-type systems, both of these venues underwent magmatism around the time of metal-logenesis (John et al., 2003; Ressel and Henry, 2006). We constructed additional models of Carlin (not presented here) that accounted for transient cooling of magmatic bodies by including spatial variations in basal heat flow, which decayed over a period of 200 k.y. (Banerjee et al., 2007). These hydrothermal models produced nearly identical fluid migration patterns and temperatures to those presented for sections C–C′ and D–D′ herein. When magmatism was represented, Banerjee et al. (2007) found that fault plane permeability had to be reduced by two orders of magnitude (from 10⁻¹¹ to 10⁻¹³ m²) to match fluid inclusion temperatures in the deposits (Hofstra and Cline, 2000). Since both values of fault permeability are realistic, it is difficult to constrain, using mathematical modeling alone, the role of magmatism during Carlin-type metallogenesis. However, given that many modern geothermal systems across the Great Basin lack a clear magmatic source, we feel that the incorporation of an interconnected permeable fault system extending to a depth of 5–6 km is a necessary condition needed to match modern thermal patterns and isotopic and geochemical observations.

It is important to realize that such high-permeability fault-controlled convective flow systems are localized in both space and time, reflecting the interplay between the inherited architecture of the Great Basin and episodic pulses of extension with and without magmatism.

This equation can represent flow induced by fluid-density gradients as well as imposed water-table gradients along the top boundary. The equation can represent changes in crustal flow patterns associated with transient changes in temperature but not changing head patterns. This implies that the flow system adjusts itself quickly to changes in fluid pressures relative to thermal processes. This is supported by analysis of hydraulic diffusivities that are more than an order of magnitude higher than thermal diffusivities.

Details of the fluid-rock isotope exchange reaction rate term (R_rk) in equation A7 are presented in the main text (see equations 1–4).

The value of the isotopic standard for oxygen (R_std) is based on standard mean ocean water composition (SMOW) and has a value of 0.0020052.

The components of the hydrodynamic dispersion tensor are calculated using equation A8. We incorporated the kinetics of quartz precipitation following the approach of Rimstidt and Barnes (1980) (see equation 7). We assumed equilibrium conditions to assign initial temperature dependent silica concentrations at the start of each model run. Details of how the silica dissolution and/or precipitation rate term (R_Qtz) are presented in the main text (see equations 5–9).

Our characterization of silica transport and deposition was also simple. We did not consider the kinetic effects of amorphous silica precipitation or other silica polymorphs (Dove and Rimstidt, 1994).

Using the above equation, Sanford (1997) and Bethke and Johnson (2002) showed how residence times can be strongly affected by hydrodynamic mixing of relatively young meteoric water within permeable aquifer units with older, connate water situated within adjacent aquitards. We used this equation to help characterize the contribution of relatively young meteoric water within closed loop convection cells during emplacement of the Carlin trend metallogenesis.

We thank Mark Coolbaugh and David Prudic for their review of our manuscript. This research was supported by the National Science Foundation Grant NSF-EAR 0809644 to Mark Person and Albert Hofstra, and by the U.S. Geological Survey Minerals Program, which funded the Metallogeny of the Great Basin project.