An active tectonic field for CO₂ storage management: the Hontomín onshore case study (Spain)

One of the concerns of underground CO₂ onshore storage is the triggering of induced seismicity and fault reactivation by the pore pressure increasing. Hence, a comprehensive analysis of the tectonic parameters involved in the storage rock formation is mandatory for safety management operations. Unquestionably, active faults and seal faults depicting the storage bulk are relevant parameters to be considered. However, there is a lack of analysis of the active tectonic strain field affecting these faults during the CO₂ storage monitoring. The advantage of reconstructing the tectonic field is the possibility to determine the strain trajectories and describing the fault patterns affecting the reservoir rock. In this work, we adapt a methodology of systematic geostructural analysis to underground CO₂ storage, based on the calculation of the strain field from kinematics indicators on the fault planes (e_y and e_x for the maximum and minimum horizontal shortening, respectively). This methodology is based on a statistical analysis of individual strain tensor solutions obtained from fresh outcrops from the Triassic to the Miocene. Consequently, we have collected 447 fault data in 32 field stations located within a 20 km radius. The understanding of the fault sets' role for underground fluid circulation can also be established, helping further analysis of CO₂ leakage and seepage. We have applied this methodology to Hontomín onshore CO₂ storage facilities (central Spain). The geology of the area and the number of high-quality outcrops made this site a good candidate for studying the strain field from kinematics fault analysis. The results indicate a strike-slip tectonic regime with maximum horizontal shortening with a 160 and 50^∘ E trend for the local regime, which activates NE–SW strike-slip faults. A regional extensional tectonic field was also recognized with a N–S trend, which activates N–S extensional faults, and NNE–SSW and NNW–SSE strike-slip faults, measured in the Cretaceous limestone on top of the Hontomín facilities. Monitoring these faults within the reservoir is suggested in addition to the possibility of obtaining a focal mechanism solutions for micro-earthquakes (M<3).

1 Introduction

Industrial human activities generate CO₂ that could change the chemical balance of the atmosphere and their relationship with the geosphere. Geological carbon storage (GSC) appears as a good choice to reduce the CO₂ gas emission to the atmosphere (Christensen and Holloway, 2004), allowing industry to increase activity with a low pollution impact. There is a lot of literature about what a site must have to be a potential underground storage suitable to GSC (e.g., Chu, 2009; Orr, 2009; Goldberg et al., 2010 among others). Reservoir sealing, caprock, permeability and porosity, plus injection pressure and volume injected, are the main considerations to choose one geological subsurface formation as the CO₂ host rock. In this context, the tectonically active field is considered in two principal ways: (1) to prevent the fault activation and earthquakes triggering, with the consequence of leakage and seepage, and (2) the long-term reservoir behavior, understanding long-term as a centennial to millennial time span. Therefore, what is the long-term behavior of GSC? What do we need to monitor for safe GSC management? Winthaegen et al. (2005) suggest three subjects for monitoring: (a) the atmosphere air quality near the injection facilities, due to the CO₂ toxicity (values greater than 4 %; see Rice, 2003, and Permentier et al., 2017), (b) the overburden monitoring faults and wells, and (c) the sealing of the reservoir. The study of natural analogs for GSC is a good strategy to estimate the long-term behavior of the reservoir, considering parameters such as the injected CO₂ pressure and volume, plus the brine mixing with CO₂ (Pearce, 2006). Hence, the prediction of site performance over long timescales also requires an understanding of CO₂ behavior within the reservoir, the mechanisms of migration out of the reservoir, and the potential impacts of a leak on the near-surface environment. The assessments of such risks will rely on a combination of predictive models of CO₂ behavior, including the fluid migration and the long-term CO₂–porewater–mineralogical interactions (Pearce, 2006). Once again, the tectonically active field interacts directly with this assessment. Moreover, the fault reactivation due to the pore pressure increasing during the injection and storage also has to be considered (Röhmann et al., 2013). Despite the uplift measured by Röhmann et al. (2013) being sub-millimeter (ca. 0.021 mm) at the end of the injection processes, given the ongoing occurrence of micro-earthquakes, long-term monitoring is required. The geomechanical and geological models predict the reservoir behavior and the caprock sealing properties. The role of the faults inside these models is crucial for the tectonic long-term behavior and the reactivation of faults that could trigger earthquakes.

Concerning induced seismicity, Wilson et al. (2017) published the Hi-Quake database, with a classification of all man-made earthquakes according to the literature, in an online repository (https://inducedearthquakes.org/, last access: May 2019). This database includes 834 projects with proved induced seismicity, where two different cases with earthquakes as large as M=1.7, detected in swarms of about 9500 micro-earthquakes, are related to GSC operations. Additionally, Foulger et al. (2018) pointed out that GSC can trigger earthquakes with magnitudes less than M=2: the cases described in their work are as great as M=1.8, with the epicenter location 2 km around the facilities. McNamara (2016) described a comprehensive method and protocol for monitoring a GSC reservoir for the assessment and management of induced seismicity. The knowledge of active fault patterns and the stress or strain field could help design a monitoring network and identify those faults capable of triggering micro-earthquakes (M<2) and/or breaking the sealing for leakage (patterns of open faults for low-permeability CO₂ migration).

In this work, we propose that the description, the analysis and establishment of the tectonic strain field have to be mandatory for long-term GSC monitoring and management, implementing the fault behavior in the geomechanical models. This analysis does not increase the cost of long-term monitoring, given that it is low-cost and the results are acquired in a few months. Therefore, we propose a methodology based on the reconstruction of the strain field from the classical studies in geodynamics (Angelier, 1979, 1984; Reches, 1983, 1987). As an innovation, we introduce the strain field (SF) analysis between 20 km away from the subsurface reservoir deep geometry, under the area of influence of induced seismicity for fluid injection. The knowledge of the strain field at a local scale allows classifying the type of faulting and their role for leakage processes, whilst the regional scale explores the tectonically active faults that could affect the reservoir. The methodology is rather simple, taking measures of slickensides and striations on fault planes to establish the orientation of the maximum horizontal shortening (e_y) and the minimum horizontal shortening (e_x) for the strain tensor. The principal advantage of the SF analysis is the direct classification of all the faults involved in the geomechanical model and the prediction of the failure parameters. Moreover, a Mohr–Coulomb failure analysis was performed on the fault pattern recognized in the Cretaceous outcrop located on top of the pilot plant.

The tectonic characterization of the GSC of Hontomín was implemented in the geological model described by Le Gallo and de Dios (2018). Beyond the use of induced seismicity and potentially active faults, the scope of this method is to propose an initial analysis to manage underground storage operations. We present how the structural analysis of fault or slip data can improve the knowledge on the tectonic large-scale fault network for the potential seismic reactivation during fluid injection and the time-dependent scale of fluid stays.

2 Hontomín onshore case study

2.1 Geological description of the reservoir

The CO₂ storage site of Hontomín is enclosed in the southern section of the Mesozoic Basque–Cantabrian Basin, known as the Burgalesa Platform (Serrano and Martínez del Olmo, 1990; Tavani, 2012), within the sedimentary Bureba Basin (Fig. 1). This geological domain is located in the northern junction of the Cenozoic Duero and Ebro basins, forming an ESE-dipping monocline bounded by the Cantabrian Mountain Thrust to the north, the Ubierna Fault System (UFS) to the south and the Asturian Massif to the west (Fig. 1).

Figure 1 (a) Location map of the study area in the Iberian Peninsula, along with the geological map of the Asturian and Basque–Cantabrian areas, labeling major units and faults (modified after Quintà and Tavani, 2012); (b) geographical location of Hontomín pilot plant (red dot) within the Basque–Cantabrian Basin. This basin is tectonically controlled by the Ubierna Fault System (UFS; NW–SE-oriented) and the parallel Polientes syncline, the Duero and Ebro Tertiary basins, and Poza de la Sal evaporitic diapir. Cret: Cretaceous; F: facies.

The Meso-Cenozoic tectonic evolution of the Burgalesa Platform starts with a first rift period during Permian and Triassic times (Dallmeyer and Martínez-García, 1990; Calvet et al., 2004), followed by relative tectonic quiescence during Early and Middle Jurassic times (e.g., Aurell et al., 2002). The main rifting phase took place during Late Jurassic and Early Cretaceous times, due to the opening of the North Atlantic and the Bay of Biscay-Pyrenean rift system (García-Mondéjar et al., 1986; Le Pichon and Sibuet, 1971; Lepvrier and Martínez-García, 1990; García-Mondéjar et al., 1996; Roca et al., 2011; Tugend et al., 2014). The convergence between Iberia and Eurasia from Late Cretaceous to Miocene times triggered the inversion of previous Mesozoic extensional faults and the development of an E–W orogenic belt (Cantabrian domain to the west and Pyrenean domain to the east) formed along the northern Iberian plate margin (Muñoz, 1992; Gómez et al., 2002; Vergés et al., 2002).

The Hontomín facilities are located within the Basque–Cantabrian Basin (Fig. 1b). The geological reservoir structure is bordered by the UFS to the south and west, by the Poza de la Sal diapir and the Zamanzas Popup structure (Carola, 2014) to the north, and by the Ebro Basin to the east (Fig. 1). The structure is defined as a forced fold-related dome structure (Tavani et al., 2013; Fig. 2), formed by an extensional fault system with migration of evaporites towards the hanging wall during the Mesozoic (Soto et al., 2011). During the tectonic compressional phase, associated with the Alpine Orogeny affecting the Pyrenees, the right-lateral transpressive inversion of the basement faults was activated, along with the reactivation of transverse extensional faults (Fig. 2; Tavani et al., 2013; Alcalde et al., 2014).

Figure 2 Interpretation of a 2-D seismic reflection profile crossing the oil exploration wells (H1, H2 and H4), along with the monitoring well (Ha) and injection well (Hi) through Hontomín pilot plant (HPP). Modified from Alcalde et al. (2014). See Fig. 1 for location; black line at the red circle.

The target reservoir and seal formations consist of Lower Jurassic marine carbonates, arranged in an asymmetric dome-like structure (Fig. 2) with an overall extent of 15 km² and located at 1485 m depth (Alcalde et al., 2013, 2014; Ogaya et al., 2013). The target CO₂ injection point is a saline aquifer formed by a dolostone unit, known as “Carniolas”, and an oolitic limestone of the Sopeña Formation, both corresponding to Lias in time (Early Jurassic). The estimated porosity of the Carniolas reaches over 12 % (Ogaya et al., 2013; Le Gallo and de Dios, 2018), and it is slightly lower at the Carbonate Lias level (8.5 % in average). The reservoir levels contain saline water with more than 20 g L⁻¹ of NaCl and very low oil content. The high porosity of the lower part of the reservoir (i.e., the Carniolas level) is the result of secondary dolomitization and different fracturing events (Alcalde et al., 2014). The minimum thickness of the reservoir units is 100 m. The potential upper seal unit comprises Lias marlstones and black shales from a hemipelagic ramp (Fig. 2; Pliensbachian and Toarcian) of the “Puerto del Pozazal” and Sopeña formations.

2.2 Regional tectonic field

The tectonic context has been described from two different approaches: (1) the tectonic style of the fractures bordering the Hontomín reservoir (De Vicente et al., 2011; Tavani et al., 2011) and (2) the tectonic regional field described from earthquakes with mechanism solutions and GPS data (Herraiz et al., 2000; Stich et al., 2006; De Vicente et al., 2008).

1. The tectonic style of the Bureba Basin was described by De Vicente et al. (2011), which classified the Basque–Cantabrian Cenozoic Basin (Fig. 1a) as transpressional with a contractional horsetail splay basin. The NW–SE-oriented Ventaniella Fault (Fig. 1a) includes the UFS in the southeastward area, being active between the Permian and Triassic period and strike-slip during the Cenozoic contraction. In this tectonic configuration, the Ubierna Fault acts as a right-lateral strike-slip fault. These authors pointed out the sharp contacts between the thrusts and the strike-slip faults in this basin. Furthermore, Tavani et al. (2011) also described a complex Cenozoic tectonic context where a right-lateral tectonic style reactivated WNW-ESE-trending faults. Both the Ventaniella and the Ubierna faults acted as transpressive structures, 120 km long and 15 km wide and featuring 0.44 mm yr⁻¹ averaged tectonic strike-slip deformation between the Oligocene and the present day. The aforementioned authors described different surface segments of the UFS of right-lateral strike-slip ranging between 12 and 14 km length. The structural data collected by Tavani et al. (2011) pointed out that 60 % of data correspond to right-lateral strike-slip with a WNW–ESE trend, together with conjugate reverse faulting with a NE–SW, NW–SE and E–W trend, and left-lateral strike-slip faults N–S-oriented. They concluded that this scheme could be related to a transpressional right-lateral tectonic system with a maximum horizontal compression, S_Hmax, striking 120^∘ E. Concerning the geological evidence of recent sediments affected by tectonic movements of the UFS, Tavani et al. (2011) suggest the middle Miocene in time for this tectonic activity. However, geomorphic markets (river and valley geomorphology) could indicate tectonic activity at present times. All of these data correspond to regional or small-scale data collected to explain the Basque–Cantabrian Cenozoic transpressive basin. The advantage of the methodology proposed here to establish the tectonic local regime affecting the reservoir is the search for local-scale tectonics (with a size of 20 km) and the estimation of the depth for the non-deformation surface for strata folding in transpressional tectonics (Lisle et al., 2009).

2. Regarding the stress field from earthquake focal mechanism solutions, Herraiz et al. (2000) pointed out the regional trajectories of S_Hmax with a NNE–SSW trend, and with a NE–SW S_Hmax trend from slip-fault inversion data. Stich et al. (2006) obtained the stress field from seismic moment tensor inversion and GPS data. These authors pointed out a NW–SE Africa–Eurasia tectonic convergence at a tectonic rate of approximately 5 mm yr⁻¹. However, no focal mechanism solutions are found within the Hontomín area (20 km) and only long-range spatial correlation could be made with high uncertainty (in time, space and magnitude). The same lack of information appears in the work of De Vicente et al. (2008), with no focal mechanism solutions in the 50 km surrounding the Hontomín pilot plant (HPP). In this work, these authors classified the tectonic regime as a uniaxial extension to strike-slip with a NW–SE S_Hmax trend.

Regional data about the tectonic field inferred from different works (Herraiz et al., 2000; Stich et al., 2006; De Vicente et al., 2008, 2011; Tavani et al., 2011; Tavani, 2012) show differences for the S_Hmax. These works explain the tectonic framework for the regional scale. Nevertheless, local tectonics could determine the low permeability and the potential induced seismicity within the reservoir. In the next section, we have applied the methodology described at Sect. 3 of this paper, in order to compare the regional results from these works.

2.3 Strategy of the ENOS European project

HPP for CO₂ onshore storage is the only one in Europe recognized as a key test facility, and it is managed and conducted by CIUDEN (Fundación Ciudad de la Energía). The HPP is located within the province of Burgos (Fig. 1b), in the northern central part of Spain.

The methodology proposed in this work and its application for long-term onshore GSC management in the context of geological risk is based on the strain tensor calculation, as part of the objectives proposed in the European project ENabling Onshore CO2 Storage in Europe (ENOS). The ENOS project is an initiative of CO2GeoNet, the European Network of Excellence on the geological storage of CO₂, to support onshore storage and highlight the associated problems, such as GSC perception, safe storage operation, potential leakage management, and health and environmental safety (Gastine et al., 2017). ENOS combines a multidisciplinary European project, which focuses in onshore storage, with the demonstration of best practices through pilot-scale projects in the case of the Hontomín facilities. Moreover, this project claims to create a favorable environment for GSC onshore through public engagement, knowledge sharing and training (Gastine et al., 2017). In this context, the work-package WP1 is devoted to “ensuring safe storage operations”.

3 Methods and rationale

The lithosphere remains in a permanent state of deformation, related to plate tectonic motion. Strain and stress fields are the consequence of this deformation on the upper lithosphere, causing different fault patterns that determine sedimentary basins and geological formations. The kinematics of these faults describe the stress or strain fields, for example measuring grooves and slickensides on fault planes (see Angelier, 1979, and Reches, 1983, among others). The relevance of the tectonic field is that stress and strain determine the earthquake occurrence by the fault activity. In this work, we have performed a brittle analysis of the fault kinematics by measuring slickenfiber on fault planes (dip–dip direction and rake), in several outcrops in the surroundings of the onshore reservoir. To carry out the methodology proposed in this work, the study area was divided into a circle with four equal areas, and we searched outcrops of fresh rock to perform the fault kinematic analysis. This allows establishing a realistic tectonic very near field to be considered during the storage seismic monitoring and long-term management. Finally, we have studied the fault plane reactivation by using the Mohr–Coulomb failure criterion (Pan et al., 2016) from the fault pattern obtained in the Cretaceous limestone outcrop located on top of the HPP facilities.

3.1 Paleostrain analysis

We have applied the strain inversion technique to reconstruct the tectonic field (paleostrain evolution), affecting the Hontomín site between the Triassic, Jurassic, Cretaceous and Neogene ages (late Miocene to present times). For a further methodological explanation, see Etchecopar et al. (1981), Reches (1983) and Angelier (1990). The main assumption for the inversion technique of fault population is the self-similarity to the scale invariance for the stress or strain tensors. This means that we can calculate the whole stress or strain fields by using the slip data on fault planes and for homogeneous tectonic frameworks. The strain tensor is an ellipsoid defined by the orientation of the three principal axes and the shape of the ellipsoid (k). This method assumes that the slip vectors, obtained from the pitch of the striation on different fault planes, define a common strain tensor or a set in a homogeneous tectonic arrangement. We assume that the strain field is homogeneous in space and time, the number of faults activated is greater than five and the slip vector is parallel to the maximum shear stress (τ).

The inversion technique is based on the Bott equations (Bott, 1959). These equations show the relationship between the orientation and the shape of the stress ellipsoid:

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }\tan \left(\mathit{\theta }\right)=[n/(l\cdot m\left)\right]\cdot \left[m^{\mathrm{2}}-\left(\mathrm{1}-n^{\mathrm{2}}\right)\cdot R^{\prime }\right],\text{(2)}& {\displaystyle }& {\displaystyle }{R}^{\prime }=\left({\mathit{\sigma }}_z-{\mathit{\sigma }}_x\right)/\left({\mathit{\sigma }}_y-{\mathit{\sigma }}_x\right),\end{array}$$

where l, m and n are the direction cosines of the normal to the fault plane, θ is the pitch of the striation and R^′ is the shape of the stress ellipsoid obtained in an orthonormal coordinate system, x, y, z. In this system, σ_y is the maximum horizontal stress, σ_x is the minimum horizontal stress axis and σ_z is the vertical stress axis.

3.2 The right-dihedral model for paleostrain analysis

The right-dihedral (RD) is a semiquantitative method based on the overlapping of compressional and extensional zones by using a stereographic plot. The final plot is an interferogram figure, which usually defines the strain regime. This method is strongly robust for conjugate fault sets and with different dip values for the same tensor. The RD was originally defined by Pegoraro (1972) and Angelier and Mechler (1977), as a geometric method, adjusting the measured fault-slip data (slickensides) in agreement with theoretical models for extension and compressive fault slip. Therefore, we can constrain the regions of maximum compression and extension related to the strain regime.

3.3 The slip model for the paleostrain analysis

The slip model (SM) is based on the Navier–Coulomb fracturing criteria (Reches, 1983), taking the Anderson model solution for this study (Anderson, 1951; Simpson, 1997). The Anderson model represents the geometry of the fault plane as monoclinic, relating the quantitative parameters of the shape parameter (K^′) with the internal frictional angle for rock mechanics (φ) (De Vicente, 1988; Capote et al., 1991). Moreover, this model is valid for neoformed faults, and some considerations have to be accounted for for previous faults and weakness planes present in the rock. These considerations are related to the dip of normal and compressional faults, such as for compressional faulting dip values lower than 45^∘, reactivated as extensional faults. This model shows the relationships between the K^′, φ and the direction cosines for the striation on the fault plane (De Vicente, 1988; Capote et al., 1991):

$$\begin{array}{}\text{(3)}& K^{\prime }=e_y/e_z,\end{array}$$

where e_z is the vertical strain axis, e_y is the maximum horizontal shortening and e_x is the minimum horizontal shortening. This model assumes that there is no change in volume during the deformation and $e_y=e_x+e_z$.

For isotropic solids, principal strain axes coincide with the principal stress axes. This means that in this work, the orientation of the principal stress axis, S_Hmax, is parallel to the orientation of the principal strain axes, e_y, and, hence, the minimum stress axis, S_hmin, is parallel to the minimum strain axis, e_x. This assumption allows us to estimate the stress trajectories (S_Hmax and S_hmin) from the e_y SM results.

Resolving the equations of Anderson for different values (Anderson, 1951), we can classify the tectonic regime that activates one fault from the measurement of the fault dip, sense of dip (0–360^∘) and pitch of the slickenside, assuming that one of the principal axes (e_x, e_y or e_z) is vertical (Angelier, 1984). We can classify the tectonic regime and represent the strain tensor by using the e_y and e_x orientation.

3.4 The K^′-strain diagram

Another analysis can be achieved by using the K^′-strain diagram developed by Kaverina et al. (1996) and codified in Python code by Álvarez-Gómez (2014). These authors have developed a triangular representation based on the fault slip, where tectonic patterns can be differentiated as strike-slip and dip-slip types. This diagram is divided into seven different zones according to the type of fault: (1) pure normal, (2) pure reverse and (3) pure strike-slip; these were combined with the possibility of oblique faults ((4) reverse strike-slip and (5) strike-slip with a reverse component) and lateral faults ((6) normal strike-slip and (7) strike-slip faults with a normal component; Fig. 3). Strike-slip faults are defined by small values for pitch (p<25^∘) and dips close to vertical planes (β>75^∘). High pitch values (p>60^∘) are related to normal or reverse fault-slip vectors. Extensional faults show e_y in the vertical, whereas compressional faults show e_y in the horizontal plane.

Figure 3 (a) Kaverina original diagram to represent the tectonic regime from an earthquake focal mechanism population (see Kaverina et al., 1996, and Álvarez-Gómez, 2014). (b) K^′-strain diagram used in this work. Dotted lines represent the original Kaverina limits. Colored zones represent the type of fault. The tectonic regime is also indicated by the relationship between the strain axes and the colored legend. SS: strike-slip. The B axis is orthogonal to the P and T axes.

This method was originally performed for earthquake focal mechanism solutions by using the focal parameters, the nodal planes (dip and strike) and rake (Kaverina et al., 1996). The triangular graph is based on the equal-area representation of the T, N, or B and P axes in spherical coordinates (T tensile, N or B neutral, and P pressure axes) and the orthogonal regression between earthquake magnitudes Ms and mb for the Harvard earthquake Global CMT (Centroid Moment Tensor) catalogue in 1996. Álvarez-Gómez (2014) presented a Python-based for computing the Kaverina diagrams, and we have modified the input parameters by including the K^′ intervals for the strain field from the SM. The relationship between the original diagram of Kaverina (Fig. 3a) and the K^′-strain diagram (Fig. 3b) that we have used in this work is shown in Fig. 3. The advantage of this diagram is the fast assignation of the type of fault and the tectonic regime that determines this fault pattern and the strain axes relationship.

Table 1 summarizes the different tectonic regimes of Fig. 3b showing the relationship with the strain main axes e_y, e_x and e_z. This diagram exhibits a great advantage to classify the type of fault according to the strain tensor. Therefore, we can assume the type of fault from the fault orientation affecting geological deposits for each strain tensor obtained.

Table 1 Different tectonic regimes, K^′ values, dip values and fault type for the Kaverina modified diagram used in this work. According to the strain axes relationship, faults can be classified and the tectonic regime can be established.

SS: strike-slip; e_x=S_hmin: minimum horizontal shortening; N: normal; e_y=S_Hmax: maximum horizontal shortening; R: reverse and e_z: vertical axis.

3.5 The circular-quadrant-search (CQS) strategy for the paleostrain analysis

In this work, we propose a low-cost strategy based on a well-known methodology for determining the stress or strain tensor affecting a GSC reservoir, which will allow the long-term monitoring of geological and seismic behavior (Fig. 4). The objective is to obtain enough structural data and spatially homogeneous of faults (Figs. 4 and 5) for reconstructing the stress or strain tensor. The key point is the determination of the orientation of the e_y, e_x and K^′ to plot on a map and, therefore, to establish the tectonic regime. We have chosen quadrants of the circles with the aim to obtain a high-quality spatial distribution of point for the interpretation of the local and very near strain field. Hence, data are homogeneously distributed, instead of being only concentrated in one quadrant of the circle.

Figure 4 Methodology proposed to obtain the strain field affecting the GSC reservoir. The distances for outcrops and quadrants proposed is 20 km. The technique of right-dihedral and the K^′-strain diagram is described in the main text. The e_y and e_x represented are a model for explaining the methodology. De_y and De_x are the direction of the maximum and minimum strain, respectively. Blue box at the center is the CO₂ storage geological underground formation.

Figure 5 Geographical location of field outcrops in the eastern part of the Burgalesa Platform domain. Black lines: observed faults; red circle: 20 km radius study zone. Rose diagram are the fault orientations from the map. A total of 447 fault data were collected in 32 outcrops. Data were measured by a tectonic compass on fault planes at outcrops. The spatial distribution of the field stations is constrained by the lithology. Coordinates are in meters, UTM H30.

Pérez-López et al. (2018) carried out a first approach to the application of this methodology at Hontomín, under the objective of the ENOS project (see next section for further details). We propose a circular searching of structural field stations (Figs. 4 and 5), located within a 20 km radius. This circle was chosen, given that active faults with the capacity of triggering earthquakes of magnitudes close to M=6 exhibit a surface rupture of tens of kilometers, according to empirical models (Wells and Coppersmith, 1994). Moreover, Verdon et al. (2015) pointed out that the maximum distance of induced earthquakes for fluid injection is 20 km. Larger distances could not be related to the stress or strain regime within the reservoir, except for the case of large geological structures (folds, master faults, etc.). Microseismicity in GSC reservoir is mainly related to the operations during the injection or depletion stages and long-term storage (Verdon, 2014; Verdon et al., 2015; McNamara, 2016).

The presence of master faults (capable of triggering earthquakes of magnitude = or > than a 6 and 5 km long segment) inside the 20 km radius circle implies that the regional tectonic field determines the strain accumulation in a kilometric fault size. Furthermore, the presence of master faults could increase the occurrence of micro-earthquakes, due to the presence of secondary faults prone to triggering earthquakes by their normal seismic cycle (Scholz, 2018). It must be borne in mind that GSC onshore reservoirs used to be deep saline aquifers (e.g., Bentham and Kirby, 2005) as in the case of Hontomín (Gastine et al., 2017; Le Gallo and de Dios, 2018), which is confined in folded and fractured deep geological structures, in which local tectonics plays a key role in micro-seismicity and the possibility of CO₂ leakage.

The constraints of this strategy are related to the absence of kinematic indicators on fault planes. This absence could occur due to later overlapping geological processes as neoformed mineralization. Also, a low rigidity eludes the slickenfiber formation, and no kinematic data will be marked on the fault plane. A poor spatial distribution of the outcrops was also taken into account for constraining the strategy. The age of sediments does not represent the age of the active deformations, and, hence, the active deformation has to be analyzed by performing alternative methods (i.e., paleoseismology, archaeoseismology).

4 Results

4.1 Strain field analysis

We have collected 447 fault-slip data points on fault planes in 32 outcrops, located within a 20 km radius circle centered at the HPP (Fig. 5). The age of the outcrops ranges from Early Triassic to post-Miocene, and they are mainly located in Cretaceous limestone and dolostone (Fig. 5, Table 2). However, no Jurassic outcrops were located, and only seven stations are located on Neogene sediments, ranging between early Oligocene to middle–late Miocene. The small number of Neogene stations is due to the mechanical properties of the affected sediments, mainly poorly lithified marls and soft-detrital fluvial deposits. Despite that, all the Neogene stations exhibit high-quality data with a number of fault-slip data ranging between 7 and 8, enough for a minimum quality analysis.

Table 2 Summary of the outcrops showing the number of faults, the type of the strain tensor obtained, the De_y, S_Hmax striking and the age of the affected geological materials. N-C is the normal component for strike-slip movement.

We have labeled the outcrops with the acronym HTM followed by a number (see Fig. 5 for the geographical location and Table 2 and Fig. 6 for the fault data). The station with the highest number of faults measured is HTM17 with 105 faults on Cretaceous limestone. Conjugate fault systems can be recognized in most of the stations (HTM1, 3, 5, 7, 10, 14, 16, 21, 23, 24, 25, 29, 30 and 32; Fig. 6), although there are a few stations with only one well-defined fault set (6, 22, 32). We have to bear in mind that the recording of conjugate fault systems is more robust for the brittle analysis than recording isolated fault sets, better constraining the solution (Žalohar and Vrabec, 2008). In total, 29 of 32 stations were used (HTM24, 27, 28 with no quality data), and of these 29 stations, 21 were analyzed with the paleostrain technique. Solutions obtained here are robust to establish the paleostrain field in each outcrop as the orientation of e_y, S_Hmax (Fig. 7).

Figure 6 Stereographic representation (cyclographic plot in Schmidt net, lower hemisphere) of the fault planes measured in the field stations. “n” is the number of available data for each geostructural station. HTM24, 27 and 28 are not included due to lack of data, and HTM17 is not included due to the high number of faults.

Figure 7 Results of the paleostrain analysis obtained and classified by age. Deym: striking of the averaged of the De_y value; F: fault stereographic representation; K^′-strain diagram with dots for each fault-slip solution; RD: right-dihedral method; SM: slip method. See the Methods section for further explanation.

The results obtained from the application of the paleostrain method have been expressed in a stereogram, RD, SM and a K^′-strain diagram (Fig. 7). The K^′-strain diagram shows the fault classification as normal faults, normal with strike-slip component, pure strike-slip, strike-slip with reverse component and reverse faults (see Fig. 3). The main faults are lateral strike-slips and normal faults, followed by reverse faults, strike-slips and oblique strike-slips faults. The results of the strain regime are as follows: (1) 43 % of extensional with shear component; (2) 22 % of shear; (3) 13 % of compressive strain (Lower Cretaceous and early–middle Miocene, Table 2); (4) 13 % of pure shear and (5) 9 % of shear with compression strain field, although with the presence of five reverse faults.

In contrast, we can observe that there are solutions with a double value for the e_y, S_Hmax orientation: HTM1, 2, 10, 11, 13, 15, 19, 26, and 30. The stations HTM3 and 23 (Upper Cretaceous) show the best solution for a strike-slip strain field as a pure strike-slip regime and e_y with a 25 and 99^∘ E trend, respectively (Fig. 7).

It is easy to observe the agreement between the e_y results from the SM and the K^′-strain diagram; for instance, in the HTM2 the K^′-strain diagram indicates strike-slip faults with a reverse component for low dips ($\mathrm{0}{}^{\circ }<\mathit{\beta }<\mathrm{40}{}^{\circ }$) but also indicates strike-slip faults with a normal component for larger dips ($\mathrm{40}{}^{\circ }<\mathit{\beta }<\mathrm{90}{}^{\circ }$). However, both results are in agreement with a strain field defined by the orientation for e_y, S_Hmax with $\mathrm{150}{}^{\circ }\pm \mathrm{18}{}^{\circ }$ trend. This tectonic field affects Cretaceous carbonates and coincides with the regional tectonic field proposed by Herraiz et al. (2000), Tavani et al. (2011) and Alcalde et al. (2014).

4.2 Late Triassic outcrop paleostrain

Strain analysis from the HTM5 fault set shows e_y with a NW–SE-trending and shear regime with the extension defined by strike-slip faults (Fig. 7a). This is in agreement with the uniaxial extension described in Tavani (2012), constraining this regime with S_Hmin with a NE–SW trend.

4.3 Cretaceous outcrop paleostrain

We have divided this result in two groups: (a) outcrops within the 20 km circle from HPP and (b) the outcrop of the HTM17 (Fig. 5), which is located in the HPP facilities and described in the next section. HTM14 is the only outcrop of Early Cretaceous age, showing a compressive tectonic stage with reverse fault solutions, defined by e_y with a NE–SW trend (Fig. 7a–c). Taking into account the extensional stage related to the main rifting stage that took place in Early Cretaceous times (i.e., Carola, 2014; Tavani, 2012; Tugend et al., 2014), we interpreted these results as a modern strain field, probably related to the Cenozoic inversion stage.

Outcrops HTM2, 3, 8, 17, 19, 20, 21, 22, 23, 25, 26, 29, 31 and 32 are from the Upper Cretaceous carbonates (Fig. 7). Results are as follows: (1) a compressive strain stage featured by e_y with a NW–SE trend, similar to the stage described in Tavani (2012), and (2) a normal strain stage with e_y striking both E–W and NE–SW (Fig. 7, HTM20, 21, 31 and 32). Finally, a (3) shear stage (activated strike-slip faults) and (4) a shear with extension (strike-slip with normal component) were described as well. These two late stages are featured by e_y with NE–SW and NW–SE trends. The existence of four different strain fields is determined by different ages during the Cretaceous and different spatial locations in relation to the main structures, the Ubierna Fault System, Hontomín Fault, Cantabrian Thrust, Montorio folded band and the Polientes syncline (Fig. 1).

Figure 8 Fault data from the outcrop HTM17 located on top of the HPP. See Fig. 5 for the geographical location. Stereogram plot is lower hemisphere and Schmidt net.

4.4 Cretaceous outcrop HTM17 on the Hontomín pilot plant

This outcrop is located on top of the geological reservoir, in a quarry of Upper Cretaceous limestones. The main advantage of this outcrop is the good development of striation and carbonate microfibers which yields high-quality data. Overall, 105 fault-slip data were measured, with the main orientation striking 75,−50 ^∘E and a conjugate set with a 130^∘ E (±10^∘) trend (Fig. 8). The result of the strain inversion technique shows an extensional field featuring an e_y trajectory striking 107^∘ E (±24^∘) related to an extensional strain field (see the K^′-strain diagram in Fig. 8). Most of the faults are extensional faults NE–SW and NW–SE-oriented (Fig. 9), in agreement with the extensional RD solution. Reverse faults are oriented NNE–SSW, E–W and WNW–ESE. The advantage of this outcrop is the geographical and stratigraphic position. It is located on top of the HPP facilities in younger sediments than the reservoir rocks. Furthermore, given that the Jurassic reservoir rock and the Cretaceous upper unit are both composed of carbonates, the fault pattern measured here could be a reflex of the fracture network affecting the Jurassic storage rocks at depth (see Figs. 2 and 9).

Figure 9 Normal and reverse faults stereograms (lower hemisphere and Schmidt net), and rose diagrams measured in HTM17. Green arrows indicate the orientation of the local paleostrain field. Grey arrows indicate the orientation of the present-day regional stress field (Herraiz et al., 2000).

4.5 Cenozoic outcrop strain field

The Cenozoic tectonic inversion in the area has been widely described by different authors (e.g., Carola, 2014; Tavani, 2012; Tungend et al., 2014). This tectonic inversion is related to compressive structures, activating NE–SW and ENE–WSW thrusts with NW–SE and NNE–SSW e_y trends, respectively. The Ubierna Fault has been inverted with right-lateral transpressive kinematics during the Cenozoic (Tavani et al., 2011). The Early Oligocene outcrop (HTM13, Fig. 7c) shows a local extensional field with e_y with a NNE–SSW and 150^∘ E trend. During the lower–middle Miocene, HTM15 and HTM30 outcrops exhibit the same e_y trend, but for a compressive tectonic regime (Fig. 7d). HTM1 shows extensional tectonics with e_y oriented 50 and 130^∘ E. Summarizing, the Cenozoic inversion and tectonic compression are detected during the early to middle Miocene and the Oligocene. However, during the middle Miocene only one extensional stage was interpreted (HTM1, Fig. 7c).

The outcrops located closer to the HPP (HTM17, 31, 32, Figs. 5 and 7) show E–W faults. HTM5 is located on the Ubierna Fault, showing a NW–SE trend, whilst HTM3 shows NE–SW strike-slip.

Strain analysis suggests that the planes parallel to the S_Hmax orientation (NNW–SSE and N–S) could induce the leakage into the reservoir (Fig. 7). Moreover, a 50^∘ E S_Hmax orientation could also affect the reservoir. HPP facilities are close to the Hontomín Fault (Fig. 5), a WNW–ESE-oriented fault, although the HTM17 station shows that N–S fault planes could play an important role for seepage of fluid into the reservoir.

5 Discussion

5.1 Regional active stress tensor in HTM17 fault pattern

The active regional field proposed by Herraiz et al. (2000), Stich et al. (2006), Tavani et al. (2011) and Alcalde et al. (2014) shows e_y, S_Hmax with almost NNW–SSE and N–S trends. Namely, the work from Herraiz et al. (2000) calculates three stress tensors within the 20 km of our study area and a Quaternary stress tensor close to the area (ca. 40 km south of Hontomín). The age of the first one is Miocene and defined by σ₁ 87∕331^∘; σ₂ 01∕151^∘; σ₃ 00∕061^∘ (dip–dip sense 0–360^∘), with an R=0.06 and S_Hmax trending 151^∘ E, under an extensional tectonic regime. Two post-Miocene stress tensors are defined by (1) σ₁ 87∕299^∘; σ₂ 00∕209^∘; σ₃ 01∕119^∘ with R=0.13, S_Hmax with a 29^∘ E trend under an extensional tectonic regime and (2) σ₁ 00∕061^∘; σ₂ 86∕152^∘; σ₃ 03∕331^∘, with R=0.76, and S_Hmax 62^∘ E under a strike-slip tectonic regime. Finally, these authors calculated a Quaternary stress tensor defined by σ₁ 85∕183^∘; σ₂ 02∕273^∘; σ₃ 03∕003^∘; R=0.02 and S_Hmax with a 101^∘ E trend under an extensional tectonic regime. The regional active stress tensor defined for Pliocene–Quaternary ages is σ₁ 88∕197^∘; σ₂ 01∕355^∘; σ₃ 00∕085^∘ for 327 data with R=0.5 and S_Hmax with an N–S trend under an extensional tectonic regional regime.

We have applied the regional active stress tensor (Herraiz et al., 2000) for studying the reactivation of previous fault patterns measured in HTM17 (Figs. 8 and 9). To carry out this study, we assume that the fault plane reactivation depends on σ₁ and σ₃ and the shape of the failure envelope. Therefore, we have used the Mohr–Coulomb failure criteria for preexisting fault planes (Xu et al., 2010; Labuz and Zang, 2012), by using the MohrPlotter v3.0 code (Allmendinger et al., 2012). Moreover, to calculate the Mohr–Coulomb circle, it is necessary to know the cohesion and friction parameters of the reservoir rock. Bearing in mind that the reservoir rocks are Lower Jurassic carbonates (dolostone and oolitic limestone; Alcalde et al., 2013, 2014; Ogaya et al., 2013), we have assumed the averaged cohesion for carbonates (limestone and dolostone) at 35^∘ and the coefficient of internal friction of 0.7 (Goodman, 1989). In addition, we have assumed no cohesion with an angle of static friction of 0.7 for preexisting faults.

Figure 10 Mohr–Coulomb failure analysis for the fault-slip data measured in HTM17 under the present-day stress tensor determined by Herraiz et al. (2000). Red dots are faults reactivated, and green and orange dots are located within the stable zone. Red rose diagram shows the orientation of reactivated faults, between N–S and 60^∘ E and from 115 to 180^∘ E. Green rose diagram shows the fault orientation for faults non-reactivated under the active stress field within the area. See text for further details. The yellow data in the Mohr–Coulomb (M–C) diagrams refer to those planes close to being reactivated and potentially reactivated by increasing the pore pressure.

Figure 10 shows the main results for the Mohr analysis. The reactivated planes under the active present stress field are red dots: 52 of the original 105 fault-slip measurements at HTM17. Green and orange dots indicate faults with no tectonic strength accumulation under the present-day stress field. Reactivated fault sets are oriented between N to 60^∘ E and 115 to 180^∘ E, with N–S and NNE–SSW as the main trends (Fig. 10, red rose diagram). Under an extensional tectonic field with R=0.5, N–S are normal faults, whereas NNE–SSW and NNW–SSE trends are strike-slips faults with extensional components. According to the results shown in Fig. 10, these faults could be reactivated without a pore pressure increase. The inactive fault orientation is constrained between 60 and 115^∘ E, mainly WNW–ESE (Fig. 10, green rose diagram). Regarding the uncertainties of these fault orientations, these values can oscillate ±5^∘, according to the field error measurement (averaged error for measuring structures by a compass).

Concerning the reliability of the results, some constraints need to be explained. The Mohr–Coulomb failure criterion is an approximation that assumes that the normal stress on the fault plane is not tensile. Furthermore, the increases in pore pressure in the reservoir rock reduces the normal stress on the plane of failure, and the interval of fault reactivation could be higher. This effect was not considered in the previous analysis since the calculation of the critical pore pressure is beyond the purpose of this work. Nevertheless, the MohrPlotter software (Allmendinger et al., 2012) allows estimating the increase of pore pressure to the critical value under some conditions.

Finally, we have applied the slip model and right-dihedral to the reactivated fault-slip data from the HTM17 outcrop (Fig. 11), by including the rake estimated from the active regional stress tensor determined by Herraiz et al. (2000). At a glance, faults oriented between 10^∘ E and 10^∘ W act as normal faults (4 out 52, Fig. 11a and c), faults between 10–50^∘ E and 10–50^∘ W act as extensional faults with strike-slip components (31 out 52), and NE–SW and NW–SE vertical faults act as pure strike-slips (8 out 52). The right-dihedral shows a tectonic regime of strike-slip with extensional component (see De Vicente et al., 1992), with orthorhombic symmetry and S_Hmax oriented 10^∘ W, which is in agreement with the stress tensor proposed by Herraiz et al. (2000) with ${\mathit{\sigma }}_{\mathrm{2}}=\mathrm{01}/\mathrm{355}$^∘ and σ₁ vertical. However, strain analysis in this case shows a strike-slip extensional tectonic regime instead of the extensional regime derived from the stress field. Despite this, both the Mohr–Coulomb analysis and the paleostrain analysis (SM and RD) suggest N–S normal faulting, NNE–SSW to NE–SW and NNW–SSE to NW–SE strike-slips as the active fault network affecting the reservoir. De Vicente et al. (1992) pointed out that the SM analysis is more robust applied to fault-slip data classified previously by other techniques. Here, we have used the Mohr–Coulomb failure criteria to separate active fault sets under the same strain tensor, yielding robustness to the results from SM and RD analysis.

Figure 11 (a) Stereogram and poles of fault sets (HTM17) reactivated under the present-day stress field suggested by Herraiz et al. (2000). (b) right-dihedral of the reactivated fault sets. (c) K^′-strain diagram showing the type of fault for each fault set.

We propose as complementary and future work a combined analysis consisting of the fault population analysis and the slip-tendency analysis (Morris et al., 1996), which could improve and identify those fault sets most likely to be reactivated under an active stress field. Although both analyses (fault population and slip tendency) are based on the stress tensor and the orientation of fault traces, the slip tendency also includes rock strength values obtained from the “in situ” tests.

5.2 Active faulting in the surroundings of HPP

Quaternary tectonic markers for the UFS are suggested by Tavani et al. (2011). According to the tectonic behavior of this fault as a right-lateral strike-slip and the fault segments proposed by Tavani et al. (2011), ranging between 12 and 14 km long, the question is whether this fault could trigger significant earthquakes and what the maximum associated magnitude could be. This is a relevant question given that the “natural seismicity” in the vicinity could affect the integrity of the caprock. Bearing in mind the expected long life for the reservoir, estimated in thousands of years, the potential natural earthquake that this master fault could trigger has to be estimated. In this sense, it is necessary to depict seismic scenarios related to large earthquake triggering; however, this type of analysis is beyond the focus of this work.

The incoming information that we have to manage in the area of influence (20 km) is (a) the instrumental seismicity, (b) the geometry of the fault, (c) the total surface rupture, (d) the upper crust thickness and (e) the heat flow across the lithosphere. Starting with the heat-flow value, the Hontomín wells show a value that lies between 62 and 78 mW m⁻² at 1500 m depth approximately (Fernández et al., 1998). Regarding the Moho depth in the area, these aforementioned authors obtained a value ranging between 36 and 40 km depth, while the lithosphere base ranges between 120 and 130 km depth (Torne et al., 2015). The relevance of this value is the study of the thermal weakness in the lithosphere that could nucleate earthquakes in intraplate areas (Holford et al., 2011). For these authors, the comparison between the crustal heat flow in particular zones, in contrast with the background regional value, could explain high seismicity and high rates of small-earthquake occurrence, as is the case of the New Madrid Seismic Zone (Landgraf et al., 2018). For example, in Australia heat-flow values of as much as 90 mW m⁻² are related to earthquakes sized M>5 (Holford et al., 2011).

Regarding the maximum expected earthquake into the zone, we have applied the empirical relationships obtained by Wells and Coppersmith (1994). We have used the equations for strike-slip earthquakes according to the strain field obtained in the area (pure shear), and the surface rupture segment for the Ubierna Fault System, assuming surface rupture segments between 12 and 14 km (Tavani et al., 2011). The obtained results show that the maximum expected earthquake ranges between M=6.0 and M=6.1. For these fault parameters, Wells and Coppersmith (1994) indicate a total area rupture ranging between 140 and 150 km². A surface fault trace rupture as lower as 7 km needs at least 20 km depth in order to reach a value of the fault area rupturing greater than 100 km², in line with a Moho between 36 and 40 at depth.

Regarding the instrumental earthquakes recorded into the area, the two largest earthquakes recorded correspond to magnitude M=3.4 and M=3.3, with a depth ranging between 8 and 11 km, respectively, and a felt macroseismic intensity of III (EMS98, http://www.ign.es/web/ign/portal, last access: May 2019). Both earthquakes occurred between 50 and 60 km of distance from the Hontomín pilot plant. Only five earthquakes have been recorded within the 20 km radius area of influence and with small magnitudes ranging between M=1.5 and M=2.3. The interesting data are the depth of these earthquakes, ranging between 10 and 20 km, which suggest that the seismogenic crust could reach 20 km depth.

5.3 Local tectonic field and induced seismicity

The fluid injection into a deep saline aquifer, which is used as GSC, generally increases the pore pressure. The increasing of the pore pressure migrates from the point of injection to the whole reservoir. Moreover, changes in the stress field for faults that are located below the reservoir could also trigger induced earthquakes (Verdon et al., 2014). Nevertheless, to understand this possibility and the study, the volumetric strain field spatial distribution is required (Lisle et al., 2009).

We have applied a physical model to estimate the total volume injected (room conditions) in reservoir conditions. Then we have applied McGarr's (2014) approximation of the maximum expected seismic moment for induced earthquakes. The injection of 10 kt of CO₂ in Hontomín (Gastine et al., 2017) represents an approximate injected volume of CO₂ of 5.56×10⁶ m³ (room conditions, pressure of 1 bar and temperature of 20 ^∘C). The P∕T conditions at the bottom of the wells have a maximum value close to 190 bar (Ortiz et al., 2015; Kovacs et al., 2015), although oscillating between 125 and 170 bar and with a maximum temperature close to 58 ^∘C. Kovacs et al. (2015) pointed out a pressure gradient of 0.023 MPa m⁻¹ and a thermal vertical gradient of 0.033 ^∘C m⁻¹, which would correspond to a pressure of 357 bar and a temperature of 51 ^∘C at 1550 m depth. P∕T bottom values obtained from the observational wells (Ha and Hi) by Ortiz et al. (2015) and Kovacs et al. (2015) were 170 bar and 42 ^∘C, respectively.

We have used the general law for gases $P_{\mathrm{1}}\cdot V_{\mathrm{1}}/T_{\mathrm{1}}=P_{\mathrm{2}}\cdot V_{\mathrm{2}}/T_{\mathrm{2}}$. Therefore, the total injected volume in reservoir conditions according to the parameters observed at the bottom of the wells are P₁=1.01 bar, T₁=20 ^∘C, $V_{\mathrm{1}}=\mathrm{5.56}\times {\mathrm{10}}^{\mathrm{6}}$ m³, P₂=170 bar and T₂=42 ^∘C. Hence, the total volume of injected CO₂ plus brine is 6.94×10⁴ m³.

McGarr (2014) empirically determined the maximum seismic moment related to a volume increasing by underground injection. The expression is Mo(max) (Nm) = G⋅ΔV (McGarr, 2014, Eq. 13), where G is the modulus of rigidity, and for the upper limit it is 3×10¹⁰ Pa, and ΔV is the total injected volume (we have used the total injected volume in reservoir conditions). The result is Mo(max) equal to 2.1×10¹⁵ Nm, which corresponds to a maximum seismic moment magnitude M_w (max) = 4.2, by applying the equation M_w = (Log Mo(max) − 9.05)∕1.5 from Hanks and Kanamori (1979), where Log is the logarithm to the base 10.

McGarr (2014) applied this approach for three cases: (1) wastewater injection, (2) hydraulic fracturing and (3) geothermal injection. We propose to include this approach for fluid injection related to geological storage of CO₂. We assume that the pore pressure increases from CO₂ injection in a similar way that wastewater does (originally defined by Frohlich, 2012). According to McGarr (2014), the utility of the analysis we have performed is “to predict in advance of a planned injection whether there will be induced seismicity”, and in the case of the HPP, to estimate the “total injected volume” in a small-scale injection plant.

Therefore, the earthquake magnitude to this fluid-injected volume according to the McGarr (2014) and Verdon et al. (2014) could be M>4 if there are faults with a minimum size of 4 km and oriented according to the present-day stress field (N–S extensional faults and NNE–SSW/NNW–SSE strike-slip faults; Fig. 10). In the case of HPP, there are faults below the reservoir with this potential earthquake triggering (Alcalde et al., 2014). Also according to McGarr (2014), this value does not have to be considered as an absolute physical limit but as a qualitative approximation. Alternatively, increasing the overpressure of the carbonate reservoir along with the pore pressure variations of about 0.5 MPa could trigger earthquakes as well. Stress drop related to fluid injections is also reported (Huang et al., 2016).

Le Gallo and de Dios (2018) described two main fault sets affecting the reservoir with N–S and E–W trends, respectively. According to the present-day stress tensor described by Herraiz et al. (2000) and Tavani et al. (2011), E–W fault sets accommodate horizontal shortening, which means that the permeability could be low. Besides, these faults are decoupled from the present-day stress tensor. However, N–S faults could act as normal faults and, hence, with a higher permeability. In this sense, the study of focal mechanism solutions could improve safety management, even for micro-earthquakes of magnitude less than M=3.

Moreover, the CO₂ lateral diffusion and pressure variation change during the fluid injection phase, and then the system would relax before being increased during the next injection phase. In this context, the intermittent and episodic injection of CO₂ could also trigger earthquakes by the stress field and fluid pressure variations in short time periods.

6 Conclusions

The application of the analysis for brittle deformation determines the active tectonic strain field, applied in the ongoing seismic monitoring for Geological Carbon Storage (GSC). The possibility that pore pressure variations due to fluid injection could change the stress or strain conditions in the reservoir's caprock, makes the study of the present-day tectonic field mandatory for storage safety operations. In this sense, we have to bear in mind that this kind of subsurface storage is designed for long-life expectancy (about thousands of years), and therefore relevant earthquakes could occur affecting the sealing and the seepage of CO₂, compromising the integrity of the reservoir. Hence, from our analysis we can conclude the following:

1. The study of the tectonic field allows classifying the geometry of the faults to prevent earthquake-prone-related structures and design a monitoring seismic network.

2. The influence area around the facilities of the GSC for studying the active stress or strain field could extend 20 km from the facility, adding missing information from the map scale and boreholes. This information could be used from the 3-D local fracture pattern estimation to avoid pore overpressure. Analysis of the stress drop due to fluid injection could be combined with this information to understand potential microseismicity associated with the injection operations.

3. In the case of the Hontomín pilot-plant, we have obtained two strain-active tectonic fields featuring shear deformation. These fields are defined by (a) a local tectonic strain field with e_y, S_Hmax striking 50^∘ E and (b) the regional one defined by e_y, S_Hmax with a 150^∘ E trend. In this context, strike-slip faults with N–S, NNE–SSW and NNW–SSE trends accumulate present-day tectonic deformation. Analysis of the Mohr–Coulomb failure criterion shows a potential reactivation of these fault sets.

4. N–S faults accumulate tectonic deformation, and they could act as normal faults. This means that this fault set is the preferential direction for potential fluid leakage. In addition, intersection with NNE–SSW and NNW–SSE could arrange 3-D networks for fluid mobilization and leakage.

5. The Ubierna Fault System represents a tectonically active fault array that could trigger natural earthquakes as large as M=6 (±0.1), from the empirical relationship of the total rupture segment (ranging between 12 and 14 km) and the total fault-area rupture (oscillating between 100 and 150 km²). Despite the lack of instrumental seismicity in the influence area, we cannot obviate the potential earthquake occurrence within intraplate areas due to the long-timescale expected life of the GSC. The heat-flow values and thermal crust conditions could determine the presence of intraplate earthquakes with magnitude M>5, for a long timescale (thousands of years), and the total injected fluid could trigger induced earthquakes greater than M=4.

6. The active strain field is now defined by an extensional tectonic defined by e_y with a N–S trend, activating N–S normal faults and right-lateral faults with NNW and NNE trends.

Finally, we state that the determination of the active tectonic strain field, the application of the slip-tendency analysis, the recognition and study of active faults within the area of influence (20 km), the estimation of the maximum potentially triggered natural earthquake, and the modeling of the stress change during the fluid injection and stress drop probably improve the operations for a secure storage. In the near future, earthquake scenarios will be the next step: modeling the Coulomb static stress changes due to fluid injection and the modeling of intensity maps of horizontal seismic acceleration.