Faults are important controls on hydrothermal circulation worldwide. More specifically, structural discontinuities, i.e. locations where faults interact and intersect, host many hydrothermal systems. In the Great Basin, western USA, an extensive characterization effort demonstrated that hydrothermal systems are controlled by one (or more) of eight types of structural discontinuities. Presumably, specific attributes of these structural settings control the generation and maintenance of permeability and porosity, and therefore localize hydrothermal processes. Herein, I examine representative examples of the eight structural settings that host hydrothermal systems in the Great Basin. For each setting, I use a boundary element method to model fault slip on the major faults and track the distribution of stress and strain in the surrounding crust. Results demonstrate that the largest magnitude and most localized stress and strain effects occur in the structural settings that host the largest number of hydrothermal systems; fault stepovers and fault terminations. Structural settings that are common in areas of strike-slip faulting also show localized stress and strain effects. The modelling presented provides process-based explanations for the empirical and conceptual results of regional characterization of Great Basin hydrothermal systems.

Thematic collection: This article is part of the Digitally enabled geoscience workflows: unlocking the power of our data collection available at: https://www.lyellcollection.org/topic/collections/digitally-enabled-geoscience-workflows

Resource assessments suggest that in addition to the existing c. 3800 MWe of electricity that are currently generated from conventional hydrothermal resources in the United States, yet to be discovered hydrothermal resources may constitute 10 times more electricity generation capacity (Williams et al. 2008, 2009; U.S. Department of Energy 2019). Many of the known and yet-to-be discovered geothermal systems in the Great Basin are so-called deep circulation fault-hosted systems (Hinz et al. 2018). To efficiently explore for these untapped resources, advancements are needed in our ability to predict the locations of new resources without the aid of extensive subsurface data. Mapping the regional distribution of structural settings that are known to host hydrothermal systems is one exploration tool with demonstrated exploration success (Craig et al. 2017, 2021; Faulds et al. 2017, 2019, 2021; Craig 2018).

In all tectonic environments, structural discontinuities (defined as locations where faults change orientation or geometry, interact, and/or intersect) are important controls on hydrothermal processes (Curewitz and Karson 1997; Rowland and Sibson 2004; Faulds et al. 2005, 2006, 2011; Cashman et al. 2012; Rowland and Simmons 2012; Faulds 2013; Moeck 2014; Faulds and Hinz 2015; Hinz et al. 2016; Wilson and Rowland 2016; Wallis et al. 2018). Deep circulation fault-hosted geothermal systems are notably present in the USA (Faulds and Hinz 2015), Indonesia (Hinz et al. 2021), the East Africa Rift (Hinz et al. 2018), and Turkey (Faulds et al. 2009), as well as a number of other places worldwide (Moeck 2014). Though this work is built upon empirical and quantitative analyses of deep circulation fault-hosted in the Great Basin these methods are directly applicable to analogous geothermal systems in the above-listed regions as well. In geothermal systems hosted by distributed permeability (e.g. Taupō Volcanic Zone, NZ or deep sedimentary basins) or other systems where deep circulation along geologic structures does not control hydrothermal processes (Moeck 2014), these methods may be less applicable.

The Great Basin, USA has been the subject of exhaustive classification of the structural settings of hydrothermal systems (Faulds et al. 2005, 2006, 2011, 2021; Faulds 2013; Faulds and Hinz 2015). As part of this effort, the structural setting was defined for 313 known hydrothermal systems (i.e. measured temperature greater than 35°C). The structural settings for hydrothermal systems in the Great Basin are fault stepovers (25% of classified systems), horse-tailing fault terminations (20%), fault intersections (18%), accommodation zones (7%), displacement transfer zones (4%), transtensional pull aparts (3%), fault bends (2%), and major normal fault zones (1%) (Faulds and Hinz 2015; Fig. 1). In many cases an individual hydrothermal system is associated with more than one of the above settings (Faulds and Hinz 2015). Of these settings, I consider all but major normal fault zones structural discontinuities because in each of the other settings the fault(s) change geometry throughout the structure and/or multiple faults intersect/interact. Major normal fault zone can therefore be considered a control case, a simple planar fault to which the other settings can be compared. Within these settings, the areas of the most significant changes (i.e. largest strike change, loci of densely spaced fault splays, fault intersections) tend to host hydrothermal systems (Fig. 1).

Kinematic and dynamic processes associated with fault slip have been invoked to explain the locations of hydrothermal upwelling relative to structural settings (Curewitz and Karson 1997; Siler et al. 2018). When slip occurs on a fault, stress is relieved along the ruptured fault segment(s) and transferred to the surrounding crust. The surrounding crust also undergoes strain (elastic and inelastic shortening and dilatation) as a result of the slip. When faults slip or dilate they may generate brecciation or become propped open by rough wall geometries or partial mineralization (Sheldon and Micklethwaite 2007). Siler et al. (2018) demonstrated that locations of maximum stress and dilatational strain affects were consistent with locations of surficial geothermal features (e.g. fumeroles, silicification, sinter), shallow temperature anomalies, and production well fields at two Great Basin geothermal fields. Though any coseismic permeability is expected to be reduced relatively rapidly via secondary mineralization (Giger et al. 2007; Zhang et al. 2008; Micklethwaite et al. 2015), recurring slip on the major faults have resulted in continued generation and maintenance of permeability at the specific locations that undergo stress change and strain (Siler et al. 2018). Over time, and with repeated earthquake cycles, these locations continually become stressed and have stresses relieved by small slip events. As a result, denser fracture networks and/or more open fractures develop relative to the surrounding crust and are therefore particularly favourable for hosting hydrothermal circulation (Siler et al. 2018).

Herein I model slip on the faults of each of the structural setting types of Faulds and Hinz (2015; Fig. 1) using a similar dynamic modelling approach to Siler et al. (2018). Results demonstrate there are inherent differences in the kinematic and dynamic effects of slip for each structural setting type. Results also indicate the relative geothermal prospectivity of each setting and define the locations where changes in stress and strain, as a result of slip, may cause long-term permeability generation and maintenance. It should be noted that the specific geometries of the faults that comprise a structural discontinuity (i.e. their strike and dip relative to regional stresses and any variations in strike or dip) may have a significant effect on the stress changes and strain that occur when the faults slip. The different structural geometries presented herein are merely representative of each structural discontinuity and the expected stress changes and strain that occur locally. It should not be assumed, however, that the models shown below are representative of all permutations of fault geometry and stress conditions.

I use a boundary element method for calculating three-dimensional (3D) displacement on triangular elements called Cut and Displace (Davis et al. 2017, 2019; Davis 2020). The modelling applies an input stress to the input faults, calculates (1) the resultant slip on the input faults, and (2) the resultant stress changes and strain that occur on fractures of a given orientation in the surrounding crust. The modelling approach allows for shear displacement on the input faults, but no opening component. Slip, strain, and stress changes are calculated for a slice at 1000 m depth. This depth represents the depth of a geothermal reservoir and is typical of reservoirs in the Great Basin (Ayling 2020). Though not universal, 1000 m depth is a reasonable simplification for our purposes here.

The geometry of the input faults for each structural setting are derived from the structural settings of Faulds and Hinz (2015; Fig. 1). Each set of fault traces is transformed into 3D assuming that normal faults dip 60° and strike-slip faults are vertical (based on simple fault mechanics (Anderson 1905)). For major normal fault zones, fault bends, horsetail fault terminations, fault stepovers, fault intersections, and accommodation zones, the six structural settings containing only normal faults, I use a normal faulting stress regime for the input stress model, for which the vertical stress (Sv) is larger than the maximum horizontal stress (SHmax), which is larger than the minimum horizontal stress (Shmin). For these structural settings, Shmin is assumed to be orthogonal to primary fault segments, azimuth 090° (east–west) for the fault geometries shown in Figures 1, 2.

Displacement transfer zones and transtensional pull aparts are structural settings that contain strike-slip faults. For these structural settings I use a strike-slip faulting stress regime for the input stress model, for which the maximum horizontal stress (SHmax), is larger than the vertical stress (Sv), which is larger than the minimum horizontal stress (Shmin). For these structural settings which contain right-lateral strike-slip faults, Shmin is oriented 60° to the strike-slip faults, azimuth 090° with respect to the fault geometries shown in Figures 1, 2.

The magnitudes of the principal stresses and pore pressure are calculated using the methods of Lund Snee and Zoback (2022). These methods assume that the crust is in a state of frictional failure equilibrium, a coefficient of friction = 0.6, and crustal density of 2600 kg m−3. 3D stress models for the strike-slip and normal faulting stress regimes are derived from the calculated stresses and pore fluid pressure (at 1000 m depth) by assuming that stress and pore pressure values at the surface are zero and linear lithostatic stress gradients with depth. All input parameters for modelling are shown in Table 1. These parameters were selected because they are representative of geological materials. The MATLAB functions used in the paper, the input data, and the results are in Siler (2023) and via ScienceBase at https://doi.org/10.5066/P9S73O5C.

Stress changes and strain that result from the modelled slip are calculated for 288 observation points with 1000 m regular spacing in a 16 km (east–west) by 18 km (north–south) grid at 1000 m depth beneath the surface. The non-terminating input faults extend 5 km beyond the edges of the observation grid, so no stress or strain effects from the edges of these faults are expected.

The stress changes are presented as Coulomb shear and normal tractions. These are calculated for fractures with the highest Coulomb shear traction increase and normal traction decrease (unclamping) for all fracture orientations within ± 1 standard deviation of the strike, dip, and rake of the input faults. In other words, the stress changes presented below are for optimally oriented fractures with similar geometry to the input faults. All strain is assumed to be elastic, and dilatation is calculated as percent volume increase. The stress effects are presented as the percent Coulomb shear traction change, and percent normal traction change from the original (pre-slip and static) stress state.

In general, along normal faults with consistent strike that are at a high angle to Shmin, dilatation, Coulomb shear stress increase, and normal stress decrease occur in the hangingwall (down-dip side) of the fault. This can be seen in the relatively straight fault segments in panels A, B, D, and E of Figures 35). At structural discontinuities, where faults are closely spaced or intersect, these simple patterns change and higher magnitude strain and stress affects are seen in some cases. Results show locations of focused dilatation (Fig. 3), focused Coulomb shear stress change (Fig. 4), and focused normal stress change (Fig. 5). In general, the highest magnitude stress and strain effects occur adjacent to and down-dip of discontinuities in the fault(s) (i.e. intersections, strike changes).

The induced dilatation for the major normal fault segment and fault bend are relatively broad and oriented along the fault strike on the downdip side, though there are subtle loci of dilatation, commonly located were the faults change strike (Fig. 3a, b). Both the major normal fault segment and fault bend undergo relatively low magnitude Coulomb shear traction (<5%) increase and normal traction decrease (>−2%) relative to the other structural settings (Figs 4a, b, 5a and b). These changes are oriented downdip and parallel to the fault traces. The magnitude of the induced dilatation (broad areas of <−0.5 × 10−3%) is relatively low compared to the horsetail and stepover, but comparable to fault intersection, displacement transfer zone, and transtensional pull apart.

The horsetail fault termination is associated with relatively localized and high magnitude dilatation (<−2.0 × 10−3%), Coulomb shear traction increase (10%), and normal traction decrease (<−8%) relative to other structural settings. These effects are focused at the northern and southern extents of the horsetail splays and on the downdip side (Figs 3c, 4c and 5c).

Similar to the fault termination, the stepover is associated with relatively localized and high magnitude dilatation (<−2.0 × 10−3%) and high magnitude Coulomb shear traction increase (>15% MPa). Dilatation and Coulomb shear traction increases are focused within and downdip of the NW–SE striking faults that bridge the north end of the stepover (Figs 3d, 4d). This same area experiences a moderate magnitude (<−4%) normal traction decrease. The decrease in normal traction occurs offset from the loci of dilatation and Coulomb shear traction increase, downdip of the intersection between the southern north–south striking normal fault and the northern most NW–SE striking fault (Fig. 5d).

Dilatation, Coulomb shear traction increases, and normal traction decreases local to the fault intersection are relatively diffuse and low-magnitude (Figs 3e, 4e and 5e), similar to the major normal fault zone and fault bend. These effects are situated along the strike of the faults, though there appears to be slight loci of dilatation and normal traction decrease in a relatively small area in the footwall (lower right) quadrant defined by the intersecting faults (Figs 3e, 5e).

The accommodation zone that I modelled initially (panel F in Figs 15) is not associated with dilatation, but rather a broad area of relatively low magnitude (<1.0 × 10−3%) shortening in the centre (Fig. 3f). Similarly, a relatively broad area of low magnitude (<1%) normal traction increase (i.e. clamping of faults) occurs in the same area (Fig. 4f). The centre of the accommodation zone is associated with a broad area of relatively low magnitude (<5% and >−5%) Coulomb shear traction increase and decrease (Fig. 5f). To further elucidate the stress changes and strain that occur at accommodation zones, two additional geometries were modelled, representing simpler end-member accommodation zones. The first is an inward-dipping accommodation zone, in which two overlapping normal faults dip towards one another (Fig. 6). The second is an outward-dipping accommodation zone, in which two overlapping normal faults dip away from one another (Fig. 7).

The inward-dipping accommodation zone is associated with moderate dilatation (<1.0 × 10−3%), minor Coulomb shear stress increase, and minor normal stress decrease between the two faults. The outward-dipping accommodation zone is associated with moderate dilatation (<1.0 × 10−3%), moderate Coulomb shear stress increase (>5%), and minor normal stress decrease downdip of the faults outside of the accommodation zone.

The displacement transfer zone is associated with a large area of moderate dilatation (>−0.5 × 10−3%). This dilatation is focused downdip of the intersection between the strike-slip fault and normal fault splays (Fig. 3g). Relatively large areas of high magnitude (>10%) Coulomb shear traction increase and moderate (<−4%) normal traction decrease also occur downdip of the intersection between the strike-slip fault and normal fault splays (Figs 4g, 5g).

The transtensional pull apart is associated with a large (c. 4+ km2) area of relatively moderate (>−0.5 × 10−3%) dilatation between the two strike-slip faults as well as south and north of the bridging normal faults. The entire area of the pull apart is also associated with moderate (>5%) Coulomb shear traction increase and low magnitude (<−5%) normal traction decrease.

Modelling results show the induced dilatation, Coulomb shear traction change, and normal traction change that occur as a result of fault slip at eight common structural settings for hydrothermal systems in the Great Basin. Modelling shows the two structural settings that most commonly host hydrothermal systems, fault stepovers (25% of Great Basin systems) and fault terminations (20%), are also clearly associated with local, relatively high magnitude dilatation, increases in Coulomb shear traction, and decreases in normal traction relative to the other structural settings (Figs 35c and d). These loci of stress and strain occur among and downdip of the bridging faults in the case of the stepover and among the horsetail fault strands in the case of the fault termination.

These coseismic stress changes and strain are short lived relative to earthquake recurrence and some of the resultant permeability is probably rapidly destroyed via secondary mineral precipitation (Giger et al. 2007; Zhang et al. 2008; Micklethwaite et al. 2015). Alternatively, faults that slip or dilate may become propped open by rough wall geometries or partial mineralization and therefore generate long-lived permeability. The dilatation, shear traction increases, and normal traction decreases suggest that fractures in these areas are more likely to be critically stressed (i.e. at or near the point of failure) and slip as result of an earthquake on the major faults. Over time and with recurring earthquake cycles, this may result in a high-density volume of faults and fractures relative to the surrounding crust, a situation that promotes local hydrothermal circulation (Siler et al. 2018). These results are largely consistent with the data and conceptual cartoons of Faulds (2013) and Faulds and Hinz (2015; Fig. 1). The modelling presented here shows the strain and stress changes are shifted downdip along the faults in both cases. In the case of the horsetail fault termination, two separate loci occur at the (northern) initiation and (southern) termination of the horsetail splays.

Stress changes and strain occurring local to the major normal fault segment, fault bend, and fault intersection scenarios are relatively low magnitude and diffuse, spread out along the fault traces. This suggests that when fault slip occurs on the major faults, there are not areas where focused changes in stress and strain are expected. Since I equate changes in stress and strain to inducing faulting and fracturing, I therefore do not expect local areas of high-density faulting and fracturing at these settings. In the case of major normal fault segments and fault bends this is largely consistent with the classification of Great Basin hydrothermal systems, with only 1 and 2% of the classified systems occurring at these structural settings, respectively (Faulds 2013; Faulds and Hinz 2015). In the case of fault intersections, which constitute 18% of Great Basin systems, the modelling herein is inconsistent with the observations of fault intersections acting as important settings for hydrothermal systems. The relatively minor stress and strain effects that occur at the fault intersection are probably a function of the specific geometry that was modelled, where one of the two faults strikes nearly parallel to Shmin. A fault intersection where both faults are better oriented for slip is likely to produce higher magnitude stress and strain results.

Modelling herein suggests that for the initial accommodation zone that was modelled (Figs 1f, 2f) shortening (not dilatation) and normal traction increases (i.e. clamping or faults rather than unclamping) occur. Neither of these effects are expected to result in slip on normal faults and promote permeability. This is surprising considering that accommodation zones host 7% of Great Basin hydrothermal systems and 26% of the systems that have been developed for electricity generation (Faulds and Hinz 2015). This result is a function of the specific accommodation zone geometry modelled, which constituted six oppositely dipping normal faults (Figs 1f, 2f). The two additional geometries that were modelled represents simpler end-member accommodation zones (Figs 6, 7). These results suggest that accommodation zones (or parts of accommodation zones) that consist of normal faults that dip toward one another are more favourable for hosting dense fracture networks and hydrothermal processes than those that consist of normal faults that dip away from one another.

The two structural settings that contain strike-slip faults, displacement transfer zone and transtensional pull apart, are both associated with relatively large areas of dilatation, Coulomb shear stress increase, and normal stress decrease. For the displacement transfer zone, this area lies downdip of the intersection between the right lateral strike-slip fault and the normal faults. For the transtensional pull apart, these effects are focused between the two right lateral strike-slip faults and downdip of the bridging normal faults. Though these two settings are not particularly common throughout the Great Basin (displacement transfer zone (4%), transtensional pull apart (3%)), they are much more common in the Walker Lane region in the western Great Basin, where strike-slip faults are far more common than to the east (Faulds et al. 2005). In the Walker Lane region displacement transfer zones and transtensional pull aparts constitute 24 and 8% of the classified systems, respectively (Faulds et al. 2005, 2006, 2011; Faulds 2013; Faulds and Hinz 2015). The modelling results are largely consistent with the data and conceptual cartoons of Faulds (2013) and Faulds and Hinz (2015), though this analysis demonstrates the largest stress and strain effects were downdip of the normal faults (which is to the SE for both geometries modelled), as opposed to the centre of the displacement transfer zone and pull apart.

The modelling presented focuses on the eight most common structural settings documented to host hydrothermal activity in the Great Basin USA. Herein, I assume the input faults are oriented ideally for reactivation relative to the input stresses. I model slip on ten different structural geometries (Figs 1, 3, 6 and 7) and calculate the resultant dilatation, Coulomb shear traction, and normal traction changes on optimally oriented faults of similar geometry (within ± 1 standard deviation of the strike, dip, and rake) to the input faults for each setting. In the areas that undergo shear traction increase, normal traction decrease, or dilation faults that slip may generate brecciation or become propped open by rough wall geometries or partial mineralization (Sheldon and Micklethwaite 2007) and therefore generate permeability.

The results reveal areas in the crust local to each structural setting where fractures are predicted to undergo dilatational strain and stress changes that may bring them closer to a critically stressed state. If the stresses overcome the frictional strength of the fracture, the fracture could slip and may generate dynamic porosity and/or permeability. If the frictional strength is not overcome, the induced stresses on the fracture may bring it closer (relative to the initial stress state) to a critically stressed-state. In such a state, the fracture is more likely to conduct fluid relative to its initial state (Barton et al. 1995; Townend and Zoback 2000; Davatzes et al. 2005; Eichhubl et al. 2009). Over long timescales (i.e. multiple earthquake cycles) volumes of the crust local to the major faults that are predicted to undergo repeated stress changes, slip, and/or dilatation are expected to develop a local high density of faults and fractures. These areas are especially well suited to host hydrothermal upwelling relative to the surrounding crust (Siler et al. 2018). Figure 8 is a summary of the results for fault bend, horsetail, stepover, fault intersection, inward dipping accommodation zone, outward-dipping accommodation zone, displacement transfer zone, and transtensional pull apart.

  • Fault stepovers and horsetail fault terminations, which host 25% and 20% of the hydrothermal systems in the Great Basin (Faulds 2013; Faulds and Hinz 2015) host the most focused and large magnitude stress changes and strain of the eight structural settings analysed. At stepovers these effects are focused within the stepover, though shifted downdip relative to the constituent faults. At fault terminations the strain and stress changes are focused near the extents of the horsetail splays to the downdip side (Fig. 8).

  • Displacement transfer zones (24% of systems in the Walker Lane strike-slip faulting region) are common structural settings for hydrothermal systems where strike-slip faults are common (Faulds 2013; Faulds and Hinz 2015). Displacement transfer zones host very large areas of high-magnitude stress changes relative (and moderate dilatation) to the other structural settings examined herein. These effects are focused outboard of the strike-slip fault and downdip of the normal faults (Fig. 8).

  • Transtensional pull aparts (8% of systems in the Walker Lane strike-slip faulting region) are also common settings for hydrothermal systems where strike-slip faults occur (Faulds 2013; Faulds and Hinz 2015). The stress change and strain effects, which are focused in and around the pull apart, are moderate relative to other settings.

  • Accommodation zones host 7% of the hydrothermal systems in the Great Basin and 26% of the systems that support electricity generation. An accommodation zone containing six faults with alternating west and east dips was originally modelled. Minor shortening (not dilatation) and normal traction increase (clamping of the faults) occurred in this setting. Simpler, two-fault accommodation zones, representing end-member inward-dipping and outward-dipping geometries, were also modelled. These results suggest that accommodation zones (or parts of accommodation zones) associated with inward-dipping faults may be more favourable for hosting hydrothermal processes (Fig. 8).

  • Similar to the accommodation zone case, fault intersections host 18% of the hydrothermal systems in the Great Basin, yet the modelling does not indicate that significant stress changes or strain are expected (at least for the single intersection geometry modelled). It is possible that different fault intersection geometries may result in more localized and focused stress and strain.

  • Major normal fault zones and fault bends are not associated with significant stress or strain localization. Accordingly, these structural settings host only 1 and 2% of Great Basin hydrothermal systems. Though not requisite, areas of more structural complexity, like stepovers, horsetail fault terminations, displacement transfer zones, transtensional pull-aparts, and inward-dipping accommodations zones are better suited to host hydrothermal systems (Fig. 8).

  • It is clear from the variety of spatial distributions and magnitudes of the induced stress changes and strain, that the geometry of the input faults is the dominant control on the effects that fault slip has on the crust. The structural settings modelled herein represent ten specific combinations of fault geometries and stress and therefore are not representative of all possible permutations. Still, the stress and strain localization patterns relative to fault geometries (Fig. 8) may be useful in determining high prospectivity areas in early–stage evaluation of geothermal areas, though examination of the specific structural setting present in any area of interest would be required. Such evaluation could be conducted using the MATLAB functions associated with this work (Siler 2023) which can be adapted to model any permutations of structure and stress.

Reviews by Cary Lindsay, Ben Melosh, and two anonymous journal reviewers significantly improved manuscript. I thank William Schermerhorn for his review of the associated MATLAB functions and data.

DLS: conceptualization (lead), data curation (lead), formal analysis (lead), funding acquisition (lead), investigation (lead), methodology (lead), software (lead), writing – original draft (lead), writing – review & editing (lead)

This work was funded by the U.S. Geological Survey.

The author declares that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

The results presented in this paper were produced by a series of MATLAB functions. Those functions, the input data and the results are available on ScienceBase via https://doi.org/10.5066/P9S73O5C.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/). Published by The Geological Society of London for GSL and EAGE.