Stress field orientation controls on fault leakage at a natural CO₂ reservoir

Travertine deposits present above the St. Johns Dome natural CO₂ reservoir in Arizona, USA, document a long (>400 kyr) history of surface leakage of CO₂ from a subsurface reservoir. These deposits are concentrated along surface traces of faults, implying that there has been a structural control on the migration pathway of CO₂-rich fluids. Here, we combine slip tendency and fracture stability to analyse the geomechanical stability of the reservoir-bounding Coyote Wash Fault for three different stress fields and two interpreted fault rock types to predict areas with high leakage risks. We find that these areas coincide with the travertine deposits on the surface, indicating that high-permeability pathways as a result of critically stressed fracture networks exist in both a fault damage zone and around a fault tip. We conclude that these structural features control leakage. Importantly, we find that even without in situ stress field data, the known leakage points can be predicted using geomechanical analyses, despite the unconstrained tectonic setting. Whilst acquiring high-quality stress field data for secure subsurface CO₂ or energy storage remains critical, we shown that a first-order assessment of leakage risks during site selection can be made with limited stress field knowledge.

1 Introduction

The successful subsurface storage of fluids in sedimentary basins is key for geo-energy technologies such as Carbon Capture and Storage (CCS), cited as a cost-effective tool for climate change mitigation, or for energy storage, required to balance the intermittency of future energy systems relying on renewable sources (Alcalde et al., 2018; Matos et al., 2019; Scott et al., 2013). The integrity of such engineered subsurface storage sites is controlled by a range of geological, geochemical and geotechnical factors. One major concern is that impermeable caprock seals may be bypassed by faults and naturally occurring, or induced, fracture networks which can form preferential fluid pathways. These could provide conduits for fluid migration, potentially leading to the rapid migration of the stored fluid (e.g. CO₂, H₂, methane) to shallow aquifers or the atmosphere (IPCC, 2005; Shipton et al., 2004; Song and Zhang, 2012). Indeed, selection criteria for subsurface storage sites commonly cite the need for minimal faulting and/or low-permeability faults intersecting or bounding the storage site (Chadwick et al., 2008; IEA GHG, 2009; Miocic et al., 2016). However, within sedimentary basins, which are key targets for geological storage of fluids, faults will occur naturally close to or within a storage complex, and thus predictability of whether a fault will act as barrier to fluid flow or not is key for an accurate risk assessment.

Whether a fault zone is sealing or non-sealing is dependent on the structure and composition of the fault zone and the mechanics of the faulting (Faulkner et al., 2010). In a widely used simple conceptual model for fault zones in siliciclastic rocks, strain is localised in the fault core that is surrounded by a damage zone of secondary structural discontinuities. Fault zones can have a single high-strain core (Chester and Logan, 1986) or contain several cores (Choi et al., 2016; Faulkner et al., 2003). The damage zone and the fault core have contrasting mechanical and hydraulic properties, with the fault core often being rich in phyllosilicates which typically have low permeability. Contrastingly, open fractures in the damage zone can have a substantially higher permeability than the host rock, if not diagenetically cemented (Caine et al., 1996; Cappa, 2009; Faulkner and Rutter, 2001; Guglielmi et al., 2008). Lateral fluid migration across the fault zone is thus controlled by (1) the permeability and continuity of the fault gouge/rock within the fault core(s), which is dependent on the host rock composition, shear strain and faulting mechanism, as well as (2) the juxtaposition of strata across the fault (Yielding et al., 1997). Inversely, vertical fluid migration is governed by fracture permeability in the damage zone (Davatzes and Aydin, 2005).

A significant amount of research has focused on understanding the mechanisms and parameters that control the composition and continuity of fault gouges as well as their permeability for different fluids as they have the potential to form effective seals (Karolytė et al., 2020; Lehner and Pilaar, 1997; Lindsay et al., 1993; Miocic et al., 2019b; Vrolijk et al., 2016). The damage zone permeability is controlled by the permeability of the host rock, the presence and geometric composition of macroscale fracture networks and deformation band networks which decrease in frequency with increasing distance from the fault core, as well as burial history, cementation and in situ stresses (Mitchell and Faulkner, 2009; Shipton et al., 2002). Outcrop studies have shown flow channelling and emphasise the strong spatial and temporal heterogeneity of fault zones (Bond et al., 2017; Burnside et al., 2013; Dockrill and Shipton, 2010; Eichhubl et al., 2009; Schulz and Evans, 1998; Soden et al., 2014). If fracture networks or faults are close to failure due to tectonically induced changes in the stress conditions or changes in pore pressure, vertical fluid flow is enhanced (Barton et al., 1995; Wiprut and Zoback, 2000). This so-called fault-valve behaviour, where faults act as highly permeable pathways for fluid discharge, is particularly likely for faults that remain active while unfavourably oriented for reactivation within the prevailing stress field (Sibson, 1990). Geomechanical parameters such as slip tendency (Morris et al., 1996) or fracture stability (Handin et al., 1963; Terzaghi, 1923) can be used to assess the potential of vertical fluid flow. The latter considers pore pressure which is a critical parameter controlling reservoir integrity not only with regards to fault weakening (Hickman et al., 1995) but also with respect to the integrity of the caprock (Caillet, 1993; Sibson, 2003).

The need for improved understanding of fracture networks and the potential of fracture reactivation and/or hydromechanically fracturing of caprock due to the injection of CO₂ has been highlighted by experiences in existing industrial CO₂ storage projects. At the Sleipner storage site, fractures in thin caprock layers appear to control the size and extent of the CO₂ plume (Cavanagh et al., 2015). The storage site of In Salah, Algeria, where between 2004 and 2011 around 4×10⁶ t of CO₂ were injected into an anticlinal structure at ∼1800 m depth, has been the focus of many studies on fracture reactivation and hydraulic fracturing of caprocks as observations at the end of the injection period suggested that pressure had migrated vertically into the caprock (Bond et al., 2013; Michael et al., 2010; Rutqvist et al., 2010; Stork et al., 2015). The existing data indicate that injection pressures were too high for the low-permeability reservoir rock and hydraulically fractured the reservoir and the lower caprock, potentially also reactivating pre-existing fracture networks related to small-scale faults (White et al., 2014).

To study how vertical fluid flow along fault zones may be related to geomechanical parameters, we examine the naturally occurring CO₂ reservoir of the St. Johns Dome, located at the border of Arizona and New Mexico. At this site, migration of fluids from the subsurface reservoir to the surface is directly linked to faults which extend through the caprock (Gilfillan et al., 2011; Moore et al., 2005), with leakage having occurred for at least 420 kyr and still ongoing (Miocic et al., 2019a; Priewisch et al., 2014). We show that leakage locations are controlled by the orientation of the reservoir bounding fault with respect to the regional stress field.

2 Geological setting

The St. Johns Dome (or Springerville–St. Johns Dome) natural CO₂ reservoir has more than 4.7×10¹⁰ m³ of recoverable CO₂ and is located on the southeastern edge of the Little Colorado River basin on the Colorado Plateau (Fig. 1) near the Arizona Transition Zone between the Basin and Range province and Rio Grande rift tectonic provinces (Bashir et al., 2011; Rauzi, 1999). It is one of 16 known naturally occurring CO₂ reservoirs on the Colorado Plateau and one of the few known naturally occurring CO₂ reservoirs worldwide where fluids leak to the surface (Gilfillan et al., 2008, 2009; Miocic et al., 2016). The CO₂ reservoir lies within a broad, NW-trending anticline that is intersected by the steeply dipping NW–SE-trending Coyote Wash Fault (Fig. 2, Moore et al., 2005; Rauzi, 1999). This major fault appears to also to form the western boundary of the productive portion of the former commercially exploited St. Johns Dome CO₂ gas field. Normal displacement across the fault ranges from less than 30 m (Salado Springs) to more than 200 m at the apex of the Cedar Mesa Anticline, 25 km SE of Salado Springs (Embid, 2009). The fault is thought to be related to Paleogene Laramide compressional tectonics which led to monoclinal folding of the Phanerozoic strata and the reactivation of older basement structures such as the Coyote Wash Fault on the Colorado Plateau (Marshak et al., 2000). The normal displacement of the fault suggests an inversion of the reverse fault related to the Basin and Range province extension starting in the early Miocene and continuing in the Pliocene as evident from displacement of Pliocene basalt flows (Embid, 2009). The Permian reservoir rocks (siltstones, sandstones and limestones), which discordantly overlie Precambrian granites (Fig. 3), are relatively shallow at 400–700 m depth, and CO₂ is present in the gas state (Gilfillan et al., 2011). Anhydrite and mudstone beds within the Permian rocks divide the reservoir vertically into several producing zones, while Triassic and Cretaceous calcareous shales and mudstones act as seals (Fig. 3). The Permian strata include, from oldest to youngest, the Supai Formation, which consists of the Amos Wash Member, Big A Butte Member, Fort Apache Member, Corduroy Member, and the San Andres Limestone and Glorieta Sandstone. A detailed geological description of the Permian rocks can be found in Rauzi (1999). The current gas–water contact is at 1425 m above sea level (a.s.l.) and the reservoir is not filled to spill. The surface rocks are mainly Triassic to Quaternary sediments, Plio-Pleistocene volcanic rocks and travertine deposits (Fig. 2). To the NW, the CO₂ reservoir is bordered by the Holbrook Basin (Harris, 2002; Rauzi, 2000) and it is closely associated with the Plio-Pleistocene Springerville volcanic field which lies just to the south and southwest of the CO₂ reservoir (Crumpler et al., 1994; Sirrine, 1958). The basaltic volcanic field consists of more than 400 individual vents and related flows, with the oldest volcanic activity dating back to around 9 Ma and the youngest flows, which can be found 8 km NW of Springerville, to about 0.3 Ma (Condit et al., 1993; Condit and Connor, 1996). As the CO₂ within the reservoir is of magmatic origin (Gilfillan et al., 2008, 2009, 2011), charging of the reservoir is thought to be the result of degassing of magma underneath the volcanic field, with CO₂ migrating along fractures and faults through the basement into the reservoir (Miocic et al., 2019a).

Figure 1 Map showing the location of natural CO₂ reservoirs and major late Cenozoic igneous rock occurrences on the Colorado Plateau and adjacent areas (after Aldrich and Laughlin, 1984; Bashir et al., 2011; Moore et al., 2005). The St. Johns Dome reservoir is located on the southern edge of the plateau, next to the Springerville volcanic field.

Figure 2 Geological map of the St. Johns Dome natural CO₂ reservoir showing the present-day extent of the CO₂ reservoir, the location of the travertine deposits, orientation of the studied faults and the location of exploration and production wells used to build the subsurface model. Structural contours indicate the top of the Fort Apache Member and illustrate the faulted anticline setting. Stress field markers indicate the azimuth of maximum horizontal stress (SH_max) after the World Stress Map in red and from volcanic clusters in blue.

Figure 3 Stratigraphic column of the St. Johns–Springerville area. Note that Cretaceous and younger deposits are often thinner than shown in this figure. CO₂ accumulations occur in the Permian strata. After Rauzi (1999) and Embid (2009).

Expression and timing of fluid flow

The travertine deposits at St. Johns Dome are an expression of CO₂-charged fluids migrating from the subsurface to the surface. Travertine formation occurs when CO₂-rich fluids outgas CO₂ as they migrate upwards to shallower depths and lower pressure, resulting in CaCO₃ supersaturation and carbonate precipitation. As such, the St. Johns Dome travertine deposits cover a surface area of more than 30 km², spread out over more than 300 km² (Figs. 2 and 4), making them one of the greatest concentrations of travertine deposits in North America. Spatially, the travertine deposits are particularly concentrated in a 10 km long zone between Salado Springs and Lyman Lake (Fig. 4, Gilfillan et al., 2011; Moore et al., 2005). This area, where present-day travertine formation occurs (Priewisch et al., 2014), is bounded by the buried Coyote Wash Fault, and the distribution of the travertine deposits and active springs suggests that the local groundwater hydrology has been influenced by the Coyote Wash Fault (Embid, 2009). Analyses of surface springs, groundwater wells and CO₂ wells with respect to the CO₂ composition, water composition and noble gas concentrations have shown that samples taken along the Coyote Wash Fault trace are influenced by waters from depth that have been enriched in mantle-derived ³He and Ca (Gilfillan et al., 2014, 2011; Moore et al., 2005). Several modelling approaches emphasise the importance of the Coyote Wash Fault for CO₂ and He migration from the Supai Formation to the surface (Allis et al., 2004; Keating et al., 2014) as, in all models, migration of gas to the surface occurs only if the fault forms a permeable conduit through the caprocks. Soil-flux measurements indicate that there is no diffuse CO₂ leakage through the caprocks, suggesting instead that faults have controlled localised fluid flow (Allis et al., 2005). In addition to the occurrences along the NE tip of the Coyote Wash Fault (cluster A), travertine mounds follow the trace of the Buttes Fault, of which the subsurface extent is not well constrained, over a distance of more than 7 km (cluster B). Travertine mounds are also found NE of the present-day extent of the CO₂ reservoir, with no clear link to other structural elements (cluster C). It is notable that there are no indications for fluid migration in the southern half of the reservoir.

Figure 4 Geological map illustrating the travertine deposits of the St. Johns Dome area. Coloured squares indicate available U–Th dating locations with respective ages in thousands of years (Miocic et al., 2019a; Priewisch et al., 2014); see Fig. 2 for legend. Four clusters of travertine mounds can be identified: cluster A, which spreads from Lyman Lake to Salado Springs, is at the tip of the Coyote Wash Fault; cluster B follows the trace of the Buttes Fault; cluster C is not related to a known structural feature and is located north of the present-day CO₂ reservoir; and cluster D is located on the crest of the Cedar Mesa Anticline. Note that ages along the Buttes Fault (cluster B) generally get younger from north to south, whilst the travertines of cluster A show a wide range of ages without an obvious spatial correlation with age.

U-series dating of the travertine mounds shows that leakage of CO₂ from the reservoir to the surface has occurred for at least 420 kyr (Fig. 4, Miocic et al., 2019a; Priewisch et al., 2014). Several of the samples analysed by Miocic et al. (2019a) fall outside the dating limitations of the U–Th method (∼500 kyr), which indicates that leakage may have occurred over much longer timescales. This is not surprising given the age of the Springerville volcanic field (earliest activity ∼9 Ma) from where the magmatic CO₂ is almost certainly sourced. Individual seeps along the Buttes Fault have lifespans of up to 200 kyr and the massive travertine platform between Salado Springs and Lyman Lake has at least a similar lifetime. Volumetric calculations indicate that the subsurface reservoir is constantly or regularly recharged, as several times the volume of CO₂ stored in the current reservoir has leaked in the past (Miocic et al., 2019a). However, due to the long timeframe of leakage recorded by the travertine deposits, only a very low percentage (0.1 %–0.001 %) of the reservoir volume (1900 Mt CO₂) has leaked annually, and thus the site could still be seen as a suitable carbon storage site from a climate mitigation point of view (Miocic et al., 2019a).

These observations illustrate that fluid migration at the St. Johns Dome occurs along fault zones, and once migration pathways have been established they are spatially fixed for long periods (>100 kyr). This is in contrast to other fault-controlled fluid migration pathways on the Colorado Plateau, for which it is suggested that these stay open only episodically for a few thousand years after rapid fault movement and subsequently heal (Frery et al., 2015). Similar cyclic reopening and healing of fractures governing fault zone permeability have been recorded by travertine deposition at other active fault zones in Italy (Brogi et al., 2010). Spatially and temporally fixed migration pathways are concerning for subsurface storage sites, and thus the processes controlling vertical fault zone permeability at the St. Johns Dome are analysed herein.

3 Methods

In order to investigate the mechanisms governing the vertical fluid flow at the St. Johns Dome, a geomechanical analysis of the Coyote Wash Fault was conducted using slip tendency and fracture stability approaches. Slip tendency (T_s) is a method that allows for a fast assessment of the relative tendency of a fault surface to undergo slip under the present-state effective stress field. It is the ratio of resolved shear stress to resolved normal stress on a surface (Morris et al., 1996):

$$T_{\mathrm{s}}={\displaystyle \frac{\mathit{\tau }}{{\mathit{\sigma }}_{\mathrm{n}}}},$$

where τ is the shear stress and σ_n the effective normal stress acting on the fault. Slip is considered likely to occur on a surface if T_s≥μ_s, with μ_s being the coefficient of static friction which is generally assumed to be 0.6 (Byerlee, 1978; Moeck et al., 2009; Sibson, 2003). Thus, using this assumption, T_s predicts slip is more likely when the friction coefficient is exceeded but does not explicitly require fault properties to be input to model (simply requiring the normal and shear stresses and an assumed coefficient of friction). Fracture stability (F_s) is the increase in pore pressure that is required to reduce the effective stresses such that a fault plane is forced into shear, tensile or hybrid failure (Handin et al., 1963; Terzaghi, 1923). Therefore, failure envelopes must be explicitly considered to model F_s which requires inputting known rock properties such as cohesion and angle of internal friction and hence fault rock composition to model.

A 3-D geological model of the St. Johns Dome was built based on published geological maps (Embid, 2009; Sirrine, 1958), well data from 37 exploration and production wells available from the Arizona oil and gas conservation commission (well logs, horizon markers) and previously published reservoir horizon map and markers (Rauzi, 1999) using Move™. Between wells, a constant stratigraphic thickness was assumed, and for the fault a dip of 70^∘ was estimated, based on previous works (Embid, 2009; Rauzi, 1999) and a 3-D dip-domain construction (Fernandez et al., 2008) of the intersection of the fault trace with the 1∕3 arcsec DEM of the 3-D elevation programme of the United States Geological Survey (USGS). The modelled fault has 6635 faces constructed as triangles from 3525 vertices. Cut-off lines were created on the fault surface by extracting the dip from a 200 m wide patch of the horizon of interest on either side of the fault and projecting this along the dip direction until it intersected with the fault (Yielding and Freeman, 2016). The current gas–water contact is at 1494 m a.s.l. (Rauzi, 1999) and is assumed to be horizontal. Due to lack of pressure data, a hydrostatic pressure gradient is assumed (0.0105 MPa m⁻¹). Geomechanical analysis of the model was conducted with industry standard software (Move™ and TrapTester^®). As no outcropping fault rocks were available, the shale gouge ratio (SGR; Yielding et al., 1997) was used as a fault rock proxy. SGR was calculated from a V_shale log of well 10-29-31, which was calculated from the gamma ray log assuming a linear response (Asquith and Krygowski, 2004). As this method only applies to siliciclastic rocks, zonal V_shale values for evaporitic sequences – 70 % shale content for anhydrite and 55 % shale for carbonates were assumed, expecting rapid fault sealing for these lithologies (Pluymakers and Spiers, 2014), or low-permeability fault rocks (Michie et al., 2018) – were used. Resulting SGR values indicate a high potential of phyllosilicate-rich fault rocks (Fig. 5). To emphasise the uncertainty regarding the fault rock composition, two different fault rocks were used for F_s calculations: a low-friction fault rock (cohesion C=0.5 MPa; coefficient of internal friction μ=0.45) and high-friction fault rock (C=0.5, μ=0.6) with rock strength values from the TrapTester^® internal database. Note that, for modelling purposes, we assume a siliciclastic sequence; however, the stratigraphic sequence also contains ∼15 % carbonate and evaporitic rocks (Fig. 3) which may have locally significant influence on the fault rock strength. T_s results are presented using stereonets, as this allows the reader to visualise how changes in the stress field orientation would influence fault stability, while F_s results are presented on a Mohr circle, as this allows a direct visualisation of how much the pore pressure needs to change to force different parts of the fault into failure. It also allows the reader to see how changes in fault rock strength could change the pore pressure needed for fault failure. For stress field data, no in situ stress measurements from the gas field itself were available; however, in addition to World Stress Map data (Heidbach et al., 2016), the nearby Springerville volcanic field can be used to derive the orientation of the horizontal stresses as presented in the following.

Figure 5 (a) Allan diagram (Allan, 1989) of the Coyote Wash Fault: dashed lines represent the hangingwall; straight lines represent footwall cut-off lines. Yellow-coloured horizons are the main reservoirs (Amos Wash Member, Big A Butte Member and Fort Apache Member of the Permian Supai Formation). (b) Shale gouge ratio plotted onto the Coyote Wash Fault plane. SGR indicates that for most of the fault phyllosilicate-rich fault rocks are likely (SGR > 30).

Stress field at the St. Johns Dome

The location of the St. Johns Dome reservoir at the margin of the Colorado Plateau and within the greater Basin and Range province has significant impact on the stress field in the study area. It is clear that the regional stress field is highly variable, as shown by the available stress field data in the vicinity of the St. Johns Dome (50 km radius, Fig. 2) from the World Stress Map (Heidbach et al., 2016) combined with a regional study on volcanic vent orientation in the Springerville volcanic field (Table 1, Connor et al., 1992). Note that the maximum horizontal stress (SH_max) from Connor et al. (1992) is based on vent clusters linearly aligned with lengths of 11 to 20 km length (Fig. 2), and that Table 1 lists them as point measurements at the centre of the cluster. To the south and southeast of the CO₂ field, the SH_max is oriented NE–SW, while west of the reservoir the SH_max orientation is highly variable, ranging from NW–SE to E–W (Fig. 2). While these data points are associated with an uncertainty of at least ±15^∘, the orientation of the stress field for the St. Johns Dome faults is difficult to constrain. A normal faulting regime (vertical stress (S_v) > maximum horizontal stress (SH_max) > minimum horizontal stress (Sh_min)) is assumed as based on the World Stress Map and works by Kreemer et al. (2010) and Wong and Humphrey (1989) for this area of the Colorado Plateau. Integration of density logs (wells 10-29-31 and 11-16-30) gives a magnitude of S_v of 23 MPa km⁻¹. Minimal horizontal stress in normal faulting regimes is typically about 65 %–85 % of the vertical stress (Hillis, 2003), which gives a magnitude of Sh_min in the range of 15 to 19.5 MPa with the magnitude of SH_max set between S_v and Sh_min.

Table 1 Table listing published stress field indicators around the St. Johns Dome. They form three clusters: IDs 2–4; 5 and 6; 7 and 8. WSM is the World Stress Map.

As the reported stress field measurements appear to form three clusters (Table 1), three different stress fields were defined (Table 2): stress field A is similar to the stress fields indicated by measurements 7 and 8, with SH_max having an azimuth of 140^∘; stress field B is oriented similar to the stress field measurements 5 and 6 with a SH_max azimuth of 100^∘; and stress field C is similar to the stress fields indicated by measurements 2 to 4 with a SH_max azimuth of 50^∘. The solitary north–south stress field measurement (ID 1) was not considered further. For the ∼ NW–SE-trending Coyote Wash Fault, these stress fields also represent the most likely (A), moderately likely (B) and least likely (C) cases for fault reactivation. Geomechanical analysis was conducted under all three defined stress fields.

Table 2 Stress fields used for the geomechanical modelling.

4 Geomechanical controls on vertical fluid migration

The results of the geomechanical analysis of the Coyote Wash Fault highlight that the orientation of the stress field has a major impact on both the slip tendency (Fig. 6) and fracture stability (Fig. 7). Slip tendency indicates that for stress field A most parts of the fault are close to failure (T_s>0.5); for stress field B, the fault is only intermediately stressed ($\mathrm{0.3}<T_{\mathrm{s}}<\mathrm{0.5}$); and for stress field C, the fault is far away from failure (T_s<0.2). Similarly, only slight increases of pore pressure are needed to force the fault into failure under stress fields A and B and a low-friction fault rock (0.95 and 1.33 MPa, respectively). The pore pressure increase needed to force the fault into failure in the case of stress field C is much higher at 6.21 MPa. Note that slip tendency in both stress fields (A and B) is higher at the NW tip of the fault than in the SE section of the fault (Fig. 6), indicating that failure is more likely to occur in the NW. This is also true for the spatial distribution of fracture stability which, for stress fields A and C (most and least likely to fail) and a low-friction fault rock, is illustrated in Fig. 7.

Figure 6 Stereoplots illustrating the slip tendency for the Coyote Wash Fault (black dots) for the three different stress fields. The fault is at least partly close to failure for stress fields A and B, while stress field C results in very low slip tendencies. Particularly the NW fault tip (indicated on the stereoplots) where travertine deposits are found at the surface is close to failure.

Figure 7 Mohr diagrams illustrating the fracture stability for the Coyote Wash Fault within the three stress fields for two types of fault rocks. Black arrows indicate the increase in pore pressure needed to force the fault into failure. Number of fault poles in each plot is 3525.

The results of the geomechanical analysis show that the bounding fault of the St. Johns Dome CO₂ reservoir is intermediately to critically stressed for two of the three modelled stress fields (A and B). For the same stress fields, a low-friction fault rock within the Coyote Wash Fault zone results in fracture stabilities which range from less than 1 to 1.33 MPa. The most critically stressed areas are located at the NW tip of the fault (Salado Springs), while the SE part of the fault is relatively stable for all studied stress fields.

The F_s values of 0.95 and 1.33 MPa for low-friction fault rock translate to an additional CO₂ column of ∼110 and ∼160 m, respectively. Currently, the reservoir is not filled to spill and the 3-D geological model indicates that the reservoir interval at the NW part of the fault could retain an additional ∼150 m of CO₂ column. Whilst a fault does not need to be critically stressed to be hydraulically conductive, the additional filling of the reservoir with a third to half more CO₂ by volume could potentially lead to fault failure and vertical fluid migration along the fault. Evidence that the reservoir has held more CO₂ in the past is provided by older travertine deposits located outside the present-day extent of the subsurface reservoir (Figs. 2 and 4; Miocic et al., 2019a) and the fact that higher paleoreservoir pressures have been implied by a geochemical study (Moore et al., 2005). These higher reservoir pressures were likely enough to bring the NW part of the fault close to failure, and we suggest that the permeability of fracture networks within the critically stressed fault damage zone was therefore increased (Barton et al., 1995; Ito and Zoback, 2000; Min et al., 2004). In order to sustain the long periods of leakage recorded by the spatially stable travertine deposition, the fault must be critically stressed for similarly long periods. Indeed, volume calculations of how much CO₂ must have leaked to the surface based on the travertine deposits show that 1–2 orders of magnitude more CO₂ was lost from the reservoir than it can hold (Miocic et al., 2019a). It is suggested that the continuous influx of magmatic CO₂ degassing from beneath the Springerville volcanic field into the reservoir caused the fault to be close to being critically stressed – a reasoning also supported by this study.

The geomechanical analysis also demonstrates that a change of the fault orientation within the stress field should not be underestimated and can lead to failure along one part of a fault, while large parts of the fault are geomechanically stable. The strike direction of the Coyote Wash Fault changes from ESE–WNW in the southern part of the fault to NW–SE in the northern section, and this change in strike is enough to render the northern section critically stressed (in two of the stress fields modelled: A and B) – with leakage pathways being the result (Fig. 9). Higher paleo gas columns within the reservoir likely contributed the forcing the fault into failure at the northern fault tip. However, some sections of the fault in the SE also have relatively low fracture stability values (Fig. 8a) which translate to only tens of metres more supported gas column than the NW section. Yet, there are no indications for past or present leakage in the SE part of the St. Johns Dome. We argue that the stress field orientation in the SE is different from the stress field orientation in the NW area of the St. Johns Dome and that as a result the fault is far from failure towards its SE tip. This is supported by stress field measurements in the vicinity of the southern edge of the reservoir (Fig. 2, Table 1), which imply a NE–SW SH_max orientation.

Figure 8 Fracture stability plotted onto the Coyote Wash Fault surface for (a) stress field A and a low-friction fault rock and (b) stress field C and a low-friction fault rock. Note the differences in the colour scale. F_s is more critical in the NW part of the fault for stress field A, while F_s for stress field C is far from critical. Dashed lines represent the hangingwall, and straight lines represent the footwall cut-off lines; see Fig. 5 for stratigraphic context.

Figure 9 Block diagram illustrating the geological setting of the St. Johns Dome area. The fault-bound CO₂ reservoir has been filled with CO₂ from depth, and leakage of CO₂ from the reservoir to the surface occurs at the NW tip of the Coyote Wash Fault. The leakage is geomechanically controlled as the stress field orientation changes along the strike of the fault as well as the strike direction of the fault. At the NW tip of the fault slip, stability values of >0.5 and fracture stability of less than 1 MPa indicate a critically stressed fault. See Fig. 2 for a complete legend, $\sim \mathrm{10}\times$ vertical exaggeration.

Vertical migration of fluids through fault and fracture networks or corridors can be classified by their location in (1) the fault damage zone, (2) at the fault tip and (3) at the crest of a fold (Ogata et al., 2014). As evidenced by travertine deposits, vertical fluid migration at the St. Johns Dome occurred at all three types of fracture networks (Fig. 2), but considerably larger volumes of fluid migrated through fracture networks linked to faults, particular at their NW tips. This indicates that, at least at this site, faults are a higher risk factor for leakage than other migration pathways such as fracture networks along the anticline structure or capillary leakage though a caprock. Based on travertine volumes, the largest volumes of leakage occurred at the NW tip of the Coyote Wash Fault, in the area between Lyman Lake and Salado Springs (Figs. 2, 4; Miocic et al., 2019a). This indicates high-permeability fracture networks within the damage zone close to the fault tip as predicted by numerical models (Backers and Moeck, 2015; Zhang et al., 2008). The lack of similar leakage pathways observed at the SE tip of the fault can be attributed to the different stress field orientation. For geological storage, in general, the occurrence of large-scale leakage at fault tips is also concerning, as displacement at fault tips usually is low and as such fault tips are not seismically resolvable and may remain undetected. Thus, seismically resolved faults could be extended beyond the normally picked extent to include the fault tips. Similarly, faults with low displacement such as the Buttes Fault, for which significant fault-related leakage has been recorded but is thought to have a maximum displacement of <25 m, may not be detectable on seismic data. This highlights the need for a good structural understanding of any geological storage site to ensure that fault tips and small faults are considered and incorporated, possibly as an additional uncertainty parameter, into the geological model.

While the geomechanical analysis highlights the role of critically stressed faults for fluid migration at the St. Johns Dome, it is missing in situ stress field data from within the CO₂ reservoir. Such data are crucial for a detailed and reliable study of fracture and fault stability (e.g. Becker et al., 2019); however, there are cases where such in situ data are missing and a geomechanical analysis may be needed (Henk, 2005). In particular, during the site selection and appraisal of subsurface storage sites, a preliminary geomechanical analysis based on existing stress field data can identify potentially critically stressed faults. The lack of in situ data can be compensated by studying several plausible stress fields (as in this study) and including uncertainties into the geomechanical analysis (Ziegler and Heidbach, 2020). For the latter, uncertainties in the stress field orientation and magnitude and in the fault orientation should be included. Statistical approaches such as Bayesian or Markov chain Monte Carlo modelling can be useful to identify uncertainty thresholds and to determine the precision by which the geomechanical parameters need to be known in order to have reliable fault and fracture stability predictions (Bao et al., 2013; Chiaramonte et al., 2008; McFarland et al., 2012).

For the geomechanical prediction of permeable fracture networks, and thus leakage pathways at the St. Johns Dome, the stress field orientation is integral. The location of the natural CO₂ reservoir at the edge of the Colorado Plateau is the likely reason for the stress field orientation change, with clear changes in crustal composition and strength in the vicinity of the St. Johns Dome (Hendricks and Plescia, 1991; Qashqai et al., 2016). The study area is also located at the intersection of the NE–NNE-trending Jemez lineament, a tectonically active zone that is characterised by Paleogene–Quaternary extension and volcanism (Fig. 1), and the ESE-trending Arizona Transition Zone (Aldrich and Laughlin, 1984; Kreemer et al., 2010). Additionally, the presence of salt deposits in the Holbrook Basin north of the study area may also impact the local stress field (Neal and Colpitts, 1997; Rauzi, 2000). The complex regional setting at the St. Johns Dome and the associated uncertainties for geomechanical modelling further highlight the need for thorough site selection criteria for engineered fluid storage sites and adequate geological data to ensure that only reservoirs with well-understood structural frameworks are chosen.

5 Implications for geological storage applications

Geomechanical modelling suggests that vertical fluid migration from the reservoir to the surface at the St. Johns Dome natural CO₂ reservoir is controlled by fracture networks in the damage zone and tip of near-critically stressed faults. We propose that regular filling of the reservoir with CO₂ from mantle sources increased the pore pressure within the reservoir and further reduced the stability of near critically stressed faults, leading to the leakage of large volumes of CO₂ over the time span of several hundred thousand years. While the leakage rates at the St. Johns Dome are low enough to render the faulted site an adequate CO₂ store for climate mitigation, similar leakage rates could impede geological storage of methane or hydrogen, particularly at onshore storage sites, based on socially and operationally acceptable criteria.

For fault-bound subsurface storage sites for CO₂ or other fluids, the history of geomechanically controlled leakage at the St. Johns Dome clearly illustrates the need for a good understanding of regional and local stress fields and faults. In particular, the stress state of faults and fault-related fracture networks prior to fluid injection needs to be well understood in order to reduce the risk of vertical fluid migration through fractured caprock. We recommended to select areas where there are no significant regional stress field changes, as these complicate geomechanical predictions. Indeed, in situ stress data from wells are key for any advanced leakage risk prognosis. To further understand the leakage mechanisms at the St. Johns Dome, geomechanical modelling of the Buttes Fault, combined with an uncertainty assessment, is recommended. More detailed dating of the travertine deposits could reveal at which part of the faults (fault tip vs. fault damage zone) failure occurred first and provide insights into the time dynamics of leakage.