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

. 10 Travertine deposits present above the St. Johns Dome natural CO 2 reservoir in Arizona, USA, document a long (>400 ka) history of surface leakage of CO 2 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 2 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 15 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 2 or energy storage remains critical, we shown that a first order assessment of leakage risks during site selection can be made 20 with limited stress field knowledge.


Introduction
The successful subsurface storage of fluids in sedimentary basins is key for GeoEnergy 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., 25 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., CO2, H2, 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 30 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 35 of the faulting (Faulkner et al., 2010). In a widely used simple conceptual model for fault zones in siliciclastic rocks, strain is localized 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 40 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 (1) by 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). 45 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 macro-scale fracture networks and deformation band networks which decrease in frequency with increasing distance from the fault core, 50 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 55 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, 60 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 CO2 has been highlighted by experiences at existing industrial CO2 storage projects.
At the Sleipner storage site, fractures in thin caprock layers appear to control the size and extent of the CO2 plume (Cavanagh et al., 2015). The storage site of In Salah, Algeria, where between 2004 and 2011 around 4 million tons of CO2 were injected 65 into an anticlinal structure at ~1,800 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). 70 To study how vertical fluid flow along fault zones may be related to geomechanical parameters we examine the naturally occurring CO2 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 ka and is 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 75 the regional stress field.

Geological setting
The St. Johns Dome (or Springerville-St. Johns Dome) natural CO2 reservoir has more than 4.7 x 10 10 m³ of recoverable CO2 and is located on the southeastern edge of the Little Colorado River Basin on the Colorado Plateau ( Fig. 1) near to the Transition Zone between the Basin and Range and Rio Grande Rift tectonic provinces (Bashir et al., 2011;Rauzi, 1999). It is one of 80 sixteen known naturally occurring CO2 reservoirs on the Colorado Plateau and one of the few known naturally occurring CO2 reservoirs world-wide where fluids are leaking to the surface (Gilfillan et al., 2008(Gilfillan et al., , 2009Miocic et al., 2016). The CO2 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 CO2 gas field. Normal displacement across the fault ranges from 85 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 extension starting in the Early Miocene and continuing in the Pliocene as evident from displacement 90 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 the CO2 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 consist of the Amos Wash Member, Big A Butte Member, Fort Apache 95 Member, Corduroy Member, and the San Andres Limestone 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 and the reservoir 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 CO2 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 south-west of the CO2 reservoir 100 (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 CO2 within the reservoir is of magmatic origin (Gilfillan et al., 2008(Gilfillan et al., , 2009(Gilfillan et al., , 2011, charging of the reservoir is thought to be the result of degassing of magma underneath the volcanic field, with CO2 migrating along fractures and faults through the basement into the reservoir (Miocic et al., 2019a). 105

Expression and timing of fluid flow
The travertine deposits at St. Johns Dome are an expression of CO2-charged fluids migrating from the subsurface to the surface. 120 Travertine formation occurs when CO2-rich fluids outgas CO2 as they migrate upwards to shallower depths and lower pressure, resulting in CaCO3 supersaturation and carbonate precipitation. As such, the St. Johns Dome travertine deposits cover a surface area of more than 30 km 2 , spread out over more than 300 km² (Figs. 2 & 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 125 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 CO2 wells with respect to the CO2 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 3 He and Ca (Gilfillan et al., 2014(Gilfillan et al., , 2011Moore et 130 al., 2005). Several modelling approaches emphasise the importance of the Coyote Wash Fault for CO2 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 cap rocks. Soil-flux measurements indicate that there is no diffuse CO2 leakage through the cap rocks, suggesting instead that faults have controlled localized fluid flow . In addition to the occurrences along the NE tip of the Coyote Wash Fault (cluster A), travertine mounds follow the trace 135 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 CO2 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. U-series dating of the travertine mounds shows that leakage of CO2 from the reservoir to the surface has occurred for at least 420 ka (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 ka), which indicates that leakage may have occurred over much longer time-scales. This is not surprising given the age of the Springerville volcanic field (earliest activity ~9 Ma) from where the 150 magmatic CO2 is almost certainly sourced. Individual seeps along the Buttes Fault have lifespans of up to 200 ka 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 CO2 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 CO2) has leaked annually and 155 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 ka). 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 few thousand years after rapid fault movement and subsequently heal (Frery et al., 2015). Similar cyclic reopening and healing of fractures 160 governing fault zone permeability has 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.

Methods
In order to investigate the mechanisms governing the vertical fluid flow at the St. Johns Dome, a geomechanical analysis of 165 the Coyote Wash Fault was conducted using slip tendency and fracture stability approaches. Slip tendency (Ts) is a method that allows for a fast assessment of the tendency of a 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) where τ is the shear stress and σn the effective normal stress acting on the fault. Slip is likely to occur on a surface if Ts≥ µs 170 with µs being the coefficient of static friction which is generally assumed to be 0.6 (Byerlee, 1978;Moeck et al., 2009;Sibson, 2003). Fracture stability (Fs) is the increase in pore pressure that is required to reduce the effective stresses such that the fault plane is forced into failure (Handin et al., 1963;Terzaghi, 1923). In contrast to slip tendency, Fs takes rock properties such as cohesion and angle of internal friction into account and thus fault rock composition needs to be known. A 3D geological model of the St. Johns Dome was built based on published geological maps (Embid, 2009;Sirrine, 1958), 175 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 3D dip-domain construction (Fernandez et al., 2008) of the intersection of the fault trace with the 1/3 arc-second DEM of the 3D elevation programme of the USGS. The modelled fault has 6635 faces constructed as 180 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 intersects with the fault (Yielding and Freeman, 2016). The current gas-water-contact is at 1494 m above sea level (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 185 available, the Shale Gouge Ratio (Yielding et al., 1997) was used as a fault rock proxy. SGR was calculated from a Vshale 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 Vshale values for evaporitic sequences (70% shale content for anhydrite, 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 190 rich fault rocks (Fig. 5). To emphasise the uncertainty regarding the fault rock composition, two different fault rocks were used for Fs calculations: clay smear (cohesion C=0.5 MPa, coefficient of internal friction µ=0.45) and phyllosilicate (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. Ts results are presented using stereonets as this here the reader can 195 readily visualise how changes in the stress field orientation would influence fault stability while Fs 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 no in situ stress measurements 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 200 as presented in the following.

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 210 (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 (SHmax) 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 south-east of the CO2 field, the SHmax is oriented NE-SW while west of the reservoir the SHmax orientation is highly variable, ranging from NW-SE to E-W (Fig. 2). While these datapoints are associated with an uncertainty 215 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 (Sv) > maximum horizontal stress (SHmax) > minimum horizontal stress (Shmin)] 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 Sv 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 220 Shmin in the range of 15 to 19.5 MPa with the magnitude of SHmax set between Sv and Shmin.
As the reported stress field measurements appear to form three clusters (Tab. 1), three different stress fields were defined (Tab. 2): stress field A is similar to the stress fields indicated by measurements 7 and 8, with SHmax having an azimuth of 140°; stress field B is oriented similar to the stress field measurements 5 and 6 with a SHmax azimuth of 100°; and stress field C is similar to the stress fields indicated by measurements 2 to 4 with a SHmax azimuth of 50°. The solitary north-south stress field 225 measurement (ID=1) was not considered further. For the ~NW-SE trending Coyote Wash Fault these stress fields also represent the most-(A), moderately-(B) and least-likely (C) cases for fault reactivation. Geomechanical analysis was conducted under all three defined stress fields.

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 235 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 (Ts>0.5), for stress field B the fault is only intermediately stressed (0.3<Ts<0.5), and for stress field C the fault is far away from failure (Ts<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 clay smear fault rock (0.95 MPa and 1.33 MPa, respectively). The pore pressure increase needed to force the fault into failure in case of stress field C is much higher at 6.21 MPa. Note that slip tendency in 240 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 clay smear fault rock is illustrated in Figure 7.
The results of the geomechanical analysis show that the bounding fault of the St. Johns Dome CO2 reservoir is intermediately to critically stressed for two of the three modelled stress fields (A and B). For the same stress fields a weak fault rock within 245 the Coyote Wash Fault zone results in fracture stabilities which range from less than 1 MPa 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 Fs values of 0.95 MPa and 1.33 MPa for a clay smear fault rock translate to an additional CO2 column of ~110 m and ~160 m, respectively. Currently the reservoir is not filled-to-spill and the 3D geological model indicates that the reservoir 250 interval at the NW part of the fault could retain an additional ~150 m of CO2 column. Thus, additional filling of the reservoir with a third to half more CO2 by volume could lead to fault failure and vertical fluid migration along the fault. Evidence that the reservoir has held more CO2 in the past is provided by older travertine deposits located outside the present day extent of the subsurface reservoir (Figs. 2 & 3, Miocic et al., 2019a) and the fact that higher paleo-reservoir pressures have been implied by a geochemical study . These higher reservoir pressures were likely enough to bring the NW part of the 255 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 on how much CO2 must have leaked to the surface based on the travertine deposits show that one to two orders of magnitude more CO2 was lost from the reservoir than it can hold (Miocic et al., 2019a). It is suggested that 260 the continuous influx of magmatic CO2 degassing from beneath the Springerville Volcanic Field into the reservoir caused the fault to be close to being critically stresseda 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; 270 A & 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 10's of meters 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 275 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, Tab. 1) which imply a NE-SW SHmax orientation.
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 280 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, 8, 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 285 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 extend to include the fault tips.
Similarly, faults with low displacement such as the Buttes Fault, for which significant fault related leakage has been recorded 290 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.

310
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 CO2 reservoir. Such data is 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 is 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 315 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 320 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 CO2 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 325 trending Jemez Lineament, a tectonically active zone that is characterized 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 330 geological data to ensure that only reservoirs with well understood structural frameworks are chosen.

Implications for geological storage applications
Geomechanical modelling suggests that vertical fluid migration from the reservoir to the surface at the St. Johns Dome natural CO2 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 CO2 from mantle sources increased the pore pressure within the reservoir and further 335 reduced the stability of near critically stressed faults, leading to the leakage of large volumes of CO2 over the time-span of several 100 kas. While the leakage rates at the St. Johns Dome are low enough to render the faulted site an adequate CO2 store for climate mitigation, similar leakage rates could socially and operationally impede geological storage of methane or hydrogen, in particularly at onshore storage sites.
For fault-bound subsurface storage sites for CO2 or other fluids the history of geomechanically controlled leakage at the St. 340 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 345 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.

Author contribution
JM and SG designed the research project which was carried out by JM with help from GJ and input from SG. JM prepared the 350 manuscript with contributions from all co-authors.

Competing interests
The authors declare that they have no conflict of interest.