Episodic fluid flow in an active fault

Fault-hosted transient fluid flow and hydrothermal activity play an important role in a wide range of geologic processes, including metallogenesis (Weis et al., 2012), oil migration (Boles et al., 2004), fault mechanics, and metamorphism (Magee and Zoback, 1993). The time scale of fluid flow is important because the longer fluid flow is active, the more it changes subsurface temperatures and alters the rock matrix. Long-term changes in fluid flow and permeability in fault zones are driven by the competition between mineral precipitation (Lowell et al., 1993) and the creation of new flow paths by fault movement and the creation of connected fractures in fault damage zones (Sibson, 1981; Eichhubl et al., 2009). While field evidence of hydrothermal minerals precipitating in faults is common (Fisher et al., 2003; Eichhubl et al., 2009), the time scale at which new flow paths are opened by fault activity and closed by cementation is largely unknown (Woodcock et al., 2007). Data on changes in permeability or the duration of fluid flow in faults are rare, especially over geological time scales (10⁴–10⁷ yr) (Ingebritsen and Gleeson, 2017).

Determining the duration of fluid flow in hydrothermal systems is often difficult because not all hydrothermal systems generate mineral deposits that can be dated, and because the uncertainty in dating hydrothermal minerals may exceed the duration of short hydrothermal events (Skinner, 1997). Several studies of fossil hydrothermal systems have demonstrated that low-temperature thermochronology offers new opportunities to date fluid flow by quantifying the thermal history of rocks in or adjacent to fluid conduits (Luijendijk, 2012; Gorynski et al., 2014; Hickey et al., 2014).

Here, we present a new approach to quantify the history of fluid flow in the active Beowawe hydrothermal system in the Basin and Range province (Nevada, United States), based on a combination of low-temperature thermochronology at high spatial resolution with a new model of conductive and advective heat flow (Luijendijk, 2019) that allows testing different hypotheses on the history of hydrothermal activity. Previous work suggests that fluid flow in the Beowawe system may have been episodic. Variations in trace elements in sinter and salinities in fluid inclusions indicate changes in fluid sources over time (Leatherman, 2010). ¹⁴C ages and stable isotope data of hydrothermal sinter deposits show a change from a light to heavy δ¹⁸O signature that coincides with an independently dated earthquake along the Malpais fault at 7500 yr B.P. (Howald et al., 2015). The most likely explanation is that the earthquake led to an increase in permeability, the opening of new flow paths, and the incorporation of a new isotopically heavy fluid source that is in equilibrium with a 5-km-deep carbonate formation (Howald et al., 2015). The ¹⁴C and stable isotope record of fluid flow goes back 20 k.y. Our objective is to quantify the long-term (>20 k.y.) history of hydrothermal activity.

The Beowawe hydrothermal system is located in north-central Nevada (Fig. 1A), southwest of the town of Elko. The area was once the second-most active geyser field in the United States, after Yellowstone National Park. However, in C.E. 1960, geyser activity ceased due to geothermal development (White, 1992). The geysers and hot springs were located in the hanging wall of the Malpais fault, an active normal fault that formed when extension rotated to a northwest-southeast direction at ca. 10–6 Ma (Zoback et al., 1981). The offset of the Malpais fault is ∼230 m (Fig. 1C). Alluvial sediments along the fault document vertical offsets of 0.7 and 2 m caused by earthquakes at 18,700 and 7500 yr B.P., respectively (Wesnousky et al., 2005). Hydrothermal activity has generated a 1600-m-long sinter deposit at the footwall of the fault with a relatively constant thickness of 65 m, consisting predominantly of opal (Zoback, 1979). The sinter terrace is underlain by middle Miocene volcanic rocks, including dacite, basaltic andesite, and basalt, that overlie Ordovician sedimentary rocks (Struhsacker, 1980). Prior to 1960, the total fluid discharge was 0.018 m³ s⁻¹, of which one-third discharged through the geysers and hot springs and two-thirds contributed to lateral discharge through alluvial sediments (Olmsted and Rush, 1987). The low salinity of hydrothermal fluids and their stable isotope composition point to a meteoric origin. The system is recharged at a topographic high northwest of the geysers and is driven by both topography and thermal convection (Olmsted and Rush, 1987; Howald et al., 2015). Low ³He/⁴He ratios suggest that there is no magmatic heat source (Banerjee et al., 2011). Geothermometers and ¹³C concentration indicate that fluids were in contact with a carbonate reservoir at temperatures averaging 230 °C at a depth of ∼5 km (Rimstidt and Cole, 1983; Day, 1987; John et al., 2003).

We collected 17 surface rock samples of middle Miocene dacites along the Malpais fault in the vicinity of the Beowawe sinter terrace (Fig. 1B; Fig. DR1 in the GSA Data Repository¹; Table 1). For four samples, the crystallization age was dated using the zircon U-Pb method. The apatite (U-Th)/He thermochronometer (AHe) was used to constrain the post-formation thermal history. For more detailed information on the methods used in this study, see the Data Repository.

We used a new advective and conductive heat flow model called Beo (Luijendijk, 2019; https://github.com/ElcoLuijendijk/beo/), which is based on the finite-element model code escript (Gross et al., 2007a, 2007b), to simulate the thermal effects of upward fluid flux in a single fault zone and lateral flow into a connected alluvial aquifer (Fig. 1D). The modeled temperature history was combined with a helium diffusion model (Meesters and Dunai, 2002a, 2002b) to calculate AHe ages. AHe ages equal crystallization age if the apatite grains have remained at low temperatures, but will tend toward zero if temperatures have exceeded 55–80 °C over long (thousands of years) time scales (Reiners et al., 2005). The model simulates only the part of the hydrothermal system where upward fluid flow and discharge take place. Flow starts at 5 km depth in a carbonate formation that is the likely origin of upward flow. The flow rate is imposed and is based on a reconstruction of the present-day fluid discharge at the surface and the flux along the fault (Olmsted and Rush, 1987). The temperature at the lower model boundary at 8 km depth was fixed at 330 °C, based on measured regional geothermal gradients (Zoback, 1979; Howald et al., 2015). We first quantified the age of the latest episode of fluid flow in this system by calibrating the model to a temperature record from Chevron–U.S. Department of Energy exploratory well 85-18 (Iovenitti and Epperson, 1981) near the western edge of the sinter terrace (Fig. 1A). We then ran models of continuous and episodic hydrothermal activity and compared modeled and measured AHe ages. Because hydrothermal activity may have varied temporally and spatially along the Malpais fault, the AHe samples were grouped along three profiles (Figs. 1A and 1B; Table 1). A number of different model experiments were conducted to test the effects of the duration of advective heating and thermal recovery on modeled AHe ages, with heating durations ranging from 0.5 to 2 k.y. and thermal recovery ranging from 2 to 9.5 k.y.

Comparison of the modeled and measured borehole temperatures (Fig. 2G) constrains the age of the latest stage of the hydrothermal activity along the Malpais fault. The model results provide the best match to the observed temperatures for a fluid flow duration of 3 k.y. (Fig. 2) using flow rates that are based on reconstructions of fluid discharge along the fault (Olmsted and Rush, 1987). Longer run times would result in too much heating of the system and would also not reproduce the overturned geothermal gradient.

The large volume of hydrothermal sinter deposits indicates that the system must have been active much longer, ∼210 k.y. (Zoback, 1979). To quantify the long-term history of fluid flow, we used low-temperature thermochronology on surface samples close to the Malpais fault. The crystallization age of the volcanic host rock of the Beowawe hydrothermal system was dated using the zircon U-Pb method as 15.6 Ma (Table DR5 in the Data Repository). Samples at distances of >40 m from the Malpais fault show AHe ages that are equal to the U-Pb age (Fig. 1B), which signifies that they have not undergone heating of temperatures exceeding ∼55–80 °C since their deposition in the Miocene. In contrast, most samples that are close to the Malpais fault show a small but consistent decrease in AHe ages that is likely the result of heating by hydrothermal fluids (Fig. 1B). The AHe data suggest that over long time scales, hydrothermal activity was focused on the Malpais fault. However, while fluid flow originated in the fault, historically the surface discharge was focused on a line of geysers on top of the sinter terrace, ∼70 m north of the fault (Fig. 1A). The focus of hydrothermal activity on the Malpais fault is consistent with the location of the oldest generation of sinter deposits that have been found close to the fault (White, 1998).

We quantified the duration of hydrothermal activity by comparing a model of advective and conductive heat flow to measured AHe ages along three profiles perpendicular to the fault (Fig. 1A). We first modeled the effects of continuous hydrothermal activity. The best fit suggests that hydrothermal activity lasted for ∼70 k.y. for profiles B and C (Fig. 3A; Fig. DR9), and ∼90 k.y. for profile A (Fig. DR8), which is shorter than the ∼210 k.y. that are required to deposit the volume of sinter present in the footwall of the fault at current fluid discharge rates and silica concentrations (Zoback, 1979). The results also indicate that the high fluid temperatures in the system would result in a full reset (AHe age = 0 Ma) of the AHe thermochronometer in a 13-m-wide zone around the fault, bounded by a 5-m-wide zone with partially reset samples, for a duration of 70 k.y. For a duration of 200 k.y., the width of the zone with fully reset AHe ages would be 45 m. However, even though the sampling was relatively dense, we did not find a single fully reset apatite (Fig. 1B; Table 1). This confirms that over longer time scales, hydrothermal activity was not continuous.

Model experiments of episodic hydrothermal activity show that periodic heating followed by thermal recovery results in only partially reset AHe ages around the fault and no fully reset samples (Fig. 3C). The best fit was obtained by 35 cycles of heating of 1 k.y. each followed by 9 k.y. of thermal recovery (Figs. 3C and 3G; Fig. DR10). A model run with heating of 2 k.y. resulted in full resetting, whereas 500 yr was not sufficient to reduce AHe ages to their observed values (Fig. 3G; Fig. DR10C). This suggests that the individual fluid pulses lasted >0.5 k.y. but <2 k.y. The duration of thermal recovery is not well constrained; a series of model runs showed that any value of 6.5 k.y. or more provided a good fit (Fig. 3G; Fig. DR10C). This means that the hydrothermal system must have been at least 255 k.y. old.

Analysis of the model error (Fig. 3G; Fig. DR10) shows that the quantified duration of individual fluid pulses of ∼1 k.y. is relatively insensitive to the inherent uncertainty of AHe data; model runs with shorter or longer fluid pulses show much higher model errors (Fig. 3G). However, the uncertainty of the total duration of fluid flow is affected by the uncertainty and dispersion of AHe ages and is estimated as ± 80 k.y. The modeled effects of hydrothermal activity on AHe ages are strongly dependent on background exhumation rates. We assumed a background exhumation rate of 1 × 10⁻⁴ m yr⁻¹, which is an average of reported regional values (see Section DR1.8 in the Data Repository). Sensitivity analysis (Fig. DR7) shows that halving the exhumation rate would double the modeled duration of hydrothermal activity.

The much better fit of episodic fluid flow models to the AHe data (Fig. 3C) compared to models of continuous fluid flow means that either fluid flow in the Beowawe system ceased between the 1-k.y.-long fluid pulses or that overall fluid flow in the Beowawe system was continuous but shifted location every ∼1 k.y. The latter explanation is favored by the distribution of ¹⁴C ages of hydrothermal sinter deposits (Howald et al., 2015), which suggests relatively constant fluid discharge over the past 20 k.y.

The 1 k.y. length of fluid pulses suggested by the AHe data is shorter than the 3 k.y. age of the latest fluid pulse based on modeled borehole temperatures. The difference may be resolved if one assumes higher flow rates at the start of the last fluid pulse and decaying permeability and flow rates over time instead of constant rates. This would have increased subsurface temperatures and reduce the reconstructed age of the last fluid pulse. At shallow levels, temperatures are limited by the boiling temperature, and higher flow rates would therefore not have affected AHe ages and the reconstructed duration of fluid pulses of 1 k.y.

The reconstructed start of hydrothermal activity of ∼250 k.y. is much younger than the 10–6 Ma age of tectonic activity of the Malpais fault. One explanation for the relatively recent appearance of the Beowawe hydrothermal system may be that continuous fault movement has gradually increased permeability of the fault damage zone, and permeability and fracture connectivity reached a threshold value that was sufficient to channel fluids along the fault from ∼5 km depth 250 k.y. ago. The formation of connected pathways in the fault damage zone would have been slowed by the formation of a clay-rich alteration zone with low permeability (Cole and Ravinsky, 1984), and by the sealing of the fault zone with silica precipitated by hydrothermal fluids.

Although the stable isotope record of Howald et al. (2015) provides a clear indication of the effects of earthquakes on opening new flow paths, not all changes in fluid flow can be attributed to earthquakes. Comparison of the modeled heat flow and the temperature record of well 85-18 provides evidence for the opening of a new flow path along the Malpais fault ∼3 k.y. ago (Fig. 2), whereas the latest earthquake was dated as 7500 yr B.P. (Wesnousky et al., 2005). One explanation may be simply that there has been a more recent earthquake on the Malpais fault that has not been detected by paleoseismic investigations. Alternatively, new flow paths may have been opened by permeability changes induced by earthquakes in nearby faults and the resulting changes in stress (Caskey and Wesnousky, 2000).

Louis and Luijendijk gratefully acknowledge support by postdoctoral startup funding (grant number 3917542) by the University of Göttingen, Germany. Support to Person by the W.M. Keck foundation (Los Angeles, California; grant 989941) is gratefully acknowledged. We also thank Volker Karius for the X-ray diffraction measurements, and Jonas Kley and David Hindle for helpful comments on an early version of the manuscript. The manuscript benefited from comments by Peter Reiners and three anonymous reviewers.