Basement structure of the Hontomín CO 2 storage site (Spain) determined by integration of microgravity and 3-D seismic data

. A multidisciplinary study has been carried out in Hontomín (Spain) to determine the basement structural setting, its geometry and the geometry of the sedimentary succession of an area aimed to be the ﬁrst Spanish pilot plant for CO 2 storage. An integration of coincident 3-D seismic results, borehole data and unpublished microgravity data were used to reproduce the deep structure and topography of the basement and to quantify the thickness of the Triassic Keuper evaporites. The subsurface structure is characterized by a half-graben setting ﬁlled with Keuper evaporites (up to 2000 m thick), forming an extensional forced fold. All data sets clearly identify two main fault systems, compartmen-talizing the main structural domain into three differentiated blocks. These faults have been interpreted to be reactivated normal faults that have led to the formation of the Hontomín dome.


Introduction
Gravity data have proven to be useful for modelling the basement topography in settings where it is overlaid by thick sedimentary successions (e.g. Døssing et al., 2014;Engen, et al., 2006;Chappell et al., 2008). This approach is sensitive to the effective removal of the gravimetric signature of the sedimentary cover and the regional, long wavelength anomaly contribution. Thus, basement topography modelling requires a good knowledge of the sedimentary succession and its con-figuration in depth. More importantly, the nature of the gravity methods allows a mutual benefit with the seismic method: gravimetric inversion profits from constraints provided by the seismic data by helping to minimize the number of possible solutions, and gravity can aid with the imaging in areas in which the seismic method is less effective, such as in subsalt areas (e.g. Contrucci et al., 2004;Filina, et al., 2015;Staedtler et al., 2014) or crystalline environments (e.g. Malehmir et al., 2009;Hedin et al., 2014).
The Hontomín structure was chosen among other options (Prado et al., 2008) to encompass the first technological development plant for CO 2 storage in Spain. Following the European Union regulation (European Union, 2009) and best practice recommendations from international experts (e.g. IEAGHG, 2008), a multidisciplinary study was carried out to assess Hontomín's suitability for geological storage of CO 2 . This included the implementation of detailed geophysical/petrophysical (Alcalde et al., 2013a, b;Ogaya et al., 2013Ogaya et al., , 2014Ugalde et al., 2013;Vilamajó et al., 2013;Rubio et al., 2011), geomechanical (Canal et al., 2013 and hydro/geochemical Permanyer et al., 2013) characterization methods.
Among this multidisciplinary approach, the acquisition, processing and interpretation of a 3-D seismic data set provided the first detailed image of the subsurface of Hontomín (Alcalde et al., 2014). However, seismic data have problems when imaging the subsalt configurations due to the high velocities of salt rock (Rousseau et al., 2003;Sava et al., 2004). This commonly results in the disappearance of the base of the salt or in an unclear image of its surface (Leveille et al., 2005). In this context, gravity data can help to constrain the contact between the salt and the underlying basement due to the density contrast between them. The salt typically has low density (2.16-2.25 kg × m −3 ), while the layers below usually have a higher density. This density contrast can be used to resolve the geometry of the salt body (Jacoby et al., 2005).
A microgravity survey was carried out at the Hontomín Technology Development Plant (Andrés, 2012) with the goal of providing additional structural constraints to areas poorly imaged from the 3-D seismics. Gravity information becomes especially important in this area where, unlike other places (Stadler et al., 2014), the seismic image was unable to provide enough information on the deepest part of the sedimentary succession and, more importantly, of the basement (Alcalde et al., 2014). The acquired high-resolution data have an excellent quality, providing a complementary approach for unravelling the subsurface structure of the area. The good control of the sedimentary sequence has allowed us to accurately remove the sedimentary cover gravity signature from the Bouguer anomaly map in order to focus on the basement structure. The gradients presented on the resulting map suggest the existence of important faults, some of which affect the basement, which are clearly related to those reported by Alcalde et al. (2013bAlcalde et al. ( , 2014. This supports their interpretation and represents a further step in the understanding of the geological evolution of the area. The integration of geophysical data sets presented in this paper is the best contribution to the knowledge of the overall geological structure of this injection site and a good example of how multidisciplinary high-resolution approaches can shed light on the shallow and/or deep investigation of geologically complex areas.

Geological setting
The Hontomín structure is located within the doubly verging Pyrenean orogen, in the south-eastern part of the Plataforma Burgalesa, to the south of the Mesozoic Basque-Cantabrian Basin (Serrano and Martínez del Olmo, 1990;Tavani and Muñoz, 2012) (Fig. 1). This area is characterized by an ESE-dipping monocline bounded to the south by the rightlateral Ubierna Fault System (Tavani et al., 2011) and by the Sierra de Cantabria Thrust to the north. Geographically, the Plataforma Burgalesa is bounded to the south by the Duero and Ebro basins and to the north by the Sierra de Cantabria.
The Plataforma Burgalesa was affected by three major deformation stages after the Permian period. The first stage corresponds to a Permian-Triassic extension (García-Mondejar et al., 1996), which reactivated wrench-fault systems (Manspeizer, 1988), giving rise to a set of ESE-WNW and E-W trending normal faults. The second event is related to the continental break-up and opening of the Bay of Biscay (García-Mondejar et al., 1996;Ziegler, 1988). This extensional episode triggered the formation of the Basque-Cantabrian Basin and the Plataforma Burgalesa. During this deformation episode, NNE-SSW faults were formed almost perpendicularly to the previous fault systems, and the Keuper rocks provided the decoupling surface that ensured different deformation styles between the sedimentary sequences above and the basement below (Tavani and Muñoz, 2012;Tavani et al., 2013). The evaporites themselves were folded with geometries that are diagnostic of extensional forced folds (Brown, 1980;Laubscher, 1982;Tavani et al., 2011Tavani et al., , 2013. Finally, the Pyrenean orogeny provided the compressional tectonic setting in which previous faults were reactivated with inverse and lateral offsets (Quintá et al., 2012;Tavani et al., 2011Tavani et al., , 2013. Among these, the Ubierna Fault System stands out as the most prominent tectonic feature in the study area. In some places, rocks of similar ages are exposed on both fault walls despite the presence of second-order faults and folds in its southern block. This, together with the well-preserved Mesozoic extensional architecture and with macro-and meso-structural data (Tavani et al., 2011), highlights an almost exclusive strike-slip behaviour during the Cenozoic inversion stage, with a much-subordinated reverse component.
The sedimentary succession of the study area ( Fig. 2) lies over a Palaeozoic basement and includes a Mesozoic cover topped by Cenozoic sediments. The stratigraphic succession begins with the Triassic Keuper facies formed by evaporites and anhydrites which are followed by Lower Jurassic anhydrites and dolomites (Pujalte et al., 2004). Above this, a succession of Lower to Middle Jurassic pelagic and hemipelagic carbonate sediments is found. Following the stratigraphic succession, Purbeck facies of the Late Jurassic to Early Cre- taceous period formed by clays, carbonates and sandstones lay unconformably above the marine Lower to Middle Jurassic rocks. Up-sequence, a siliciclastic succession is found, comprised by the Weald facies and the Escucha and Utrillas formations, which are interpreted as fluvial successions with sandstone infilled channels alternating with flood plain sediments. The stratigraphic succession ends up with Upper Cretaceous carbonates and lacustrine and detritic Cenozoic sediments lying unconformably above the Mesozoic series (Vera, 2004).

Seismic and well data
A 3-D seismic data set was acquired in summer 2010 in Hontomín, across an area of 36 km 2 centred according to the target dome-shaped structure (Alcalde et al., 2013a) (Fig. 2). The acquired seismic data set was processed down to 1500 ms up to post-stack time migration (Alcalde et al., 2013b). Critical steps in processing were related to the existence of an unexpected velocity inversion near the surface, which decreased the general quality of the data. However, the final migrated volume provided clear images from the subsurface structures down to the anhydrite level (i.e. the bottom of the dome structure). Further details on the seismic data acquisition and processing can be found in Alcalde et al. (2013aAlcalde et al. ( , b, 2014. The seismic data were integrated with well-log data from four hydrocarbon exploration boreholes (H1 to H4) and three shallow groundwater sampling wells (GW1 to GW3) (Fig. 2). Results indicate that the structure of the Hontomín dome includes eight different sedimentary packages, from Triassic to Cenozoic and four sets of faults. The target Mesozoic dome structure has a NW-SE orientation and is bounded by two major faults, the south (S) and the east (E) faults, showing vertical offsets of up to 450 and 250 m respectively. Another two sets of faults have been singled out: they trend N-S and E-W and have been related to extensional episodes that occurred during the opening of the Bay of Biscay. The dips of the dome flanks decrease upwards suggesting a protracted and discontinuous growth. The joint interpretation of the seismic and well log data allowed for inferring three evolutionary stages for the Hontomín structure (Alcalde et al., 2104): (1) the development of N-S trending faults as a consequence of differential loading around Triassic E-W faults could have allowed the development of the lowermost Jurassic units and probably produced forced folding of the Hontomín dome.
(2) Later on, the development of E-trending faults during the opening of the Bay of Biscay is recorded by the accumulation of Purbeck deposits. Simultaneously, the migration of the Triassic salts would have further increased

Gravity data
Microgravity data were acquired between August and December 2010 by Implemental System Company. A highdensity mesh covering an area of 4 × 4 km 2 was recorded with measurements taken at every 100 m, which resulted in 1600 total measured points (Fig. 2). The measurements were carried out using a Scintrex CG-5 2008 gravity meter with a resolution of 0.001 mGal. The acquisition parameters ensured a resolution lower than 4 µGal, with over three measurements per grid point on average. Seventy-one of those gravimetric measurement points were concreted, along with the seven gravimetric stations used to calibrate the gravitymeter, delivering an accuracy of 0.5 cm. The positioning of the 1600 measured points was performed with a LEICA GPS9000 and has an accuracy of 1cm or less. The calibration of the topographic instrument was carried out with respect to a geodesic vertex located within the study area. This ensures that measurement points were within the high-quality range of the device.
The gravity data were provided without the effects of moon/sun/earth tides and instrumental drift, which were al-ready corrected by the contractor. The data were then processed using Geosoft Oasis Montaj ™ to produce the complete Bouguer anomaly (BA) map of the area. This includes carrying out the Bouguer slab, free-air and terrain corrections (Blakely, 1995). The resulting Bouguer anomaly map is shown in Fig. 3. The reduction density used was 2400 kg m −3 , although several values were tested to assess the better reduction density for performing the Bouguer slab correction. Anomaly maps obtained using densities of 2000, 2400 and 2600 kg m −3 were compared with the topography, being 2400 kg m −3 which showed less correlation with the relief.
Several procedures may be applied to Bouguer anomaly data in order to enhance features that may later on help in the interpretation. Accordingly, the following processes were carried out in the present data set: (1) separation of the long and short wavelength components, (2) calculation of the gravity derivatives (vertical and total horizontal derivatives, THD), (3) construction of two 2-D gravity models and (4) inversion of gravity data. To carry out a qualitative assessment of the source depth of the Bouguer anomaly data, the short and long wavelength number components had to be separated. This separation was performed by applying an upward continuation to the gravity field (Jacobsen, 1987) of up to 350 m. The filter was selected in order to minimize distortion and ringing effects (Fig. 4a). The regional component was then subtracted from the complete Bouguer anomaly to generate a residual anomaly map presented in Fig. 4b.
Edge enhancements are often used to aid in the potential field data interpretation (Verduzco et al., 2004). Accordingly, two sets of derivatives were performed to our data set aiming to highlight the gradient zones and accordingly the structural setting of the area. The vertical derivative ( Fig. 4c) (Blakely and Simpson, 1986) and total horizontal derivative (THD, Fig. 4d) show a coherent result between them. Both maps highlight coincident areas with maximum gradient, that is, areas where the contacts between materials of different density exist. However, the vertical derivative also enhances those that are shallow.
Finally, two 2-D gravity profiles were extracted from the Bouguer anomaly map in order to produce a model of the structure of the area along them (see location in Fig. 3): Model 1, striking NE-SW and Model 2 striking almost E-W (Fig. 5). First, these models were computed with the assistance of the information provided by wells H1 and H2 but without the input of 3-D seismic data (Andrés, 2012) (Table 1). The models were then revisited to include the variations in layer depth and fault offset observed in the seismic model (Fig. 5a). It is worth highlighting that changes introduced to the models were not significant and provide a good match with the seismic model and 2-D sections derived from the inversion procedure (Fig. 6).

Gravity inversion
Gravity data can be inverted in order to obtain the topography of a desired layer (Oldenburg, 1974). However, inversion results are more reliable when the number of variables to in-vert decrease. Inverting just the topography of one layer can provide us with excellent results if we can preliminarily establish the density of every layer and the thickness of every density interval. The method utilized in this study is based on the Parker's algorithm (Blakely, 1996;Parker, 1972). It works in the frequency domain and is integrated in GM-SYS 3-D platform. Here, we inverted the microgravity data in order to obtain the geometry of the basement. The large amount of information available in the Hontomín area has granted us the basis to build a wellconstrained model that allows us to perform a successful inversion. Well data, seismic horizons and 2-D gravity models were used to construct an initial model. The inversion process, which was then applied, consisted of isolating all the gravity contributors, i.e. different sources (sediments and regional component) and calculating their gravity responses. That was then subtracted from the Bouguer gravity anomaly (BA) map (Fig. 3) and the residual was inverted.
Accordingly, building the gravity model for the inversion consisted of the following steps: -Calculation and removal from the BA map of the gravity anomaly generated by the stratigraphic succession below 0 m ( Fig. 7a and c), assuming a constant density for each layer (Table 1).
-Analysis, calculation and removal of the long wavelength/deep signal corresponding to first-order polynomial ( Fig. 7b and d). This filter was selected in order to avoid removing the signal from the Keuper succession.
-Generation of the basement topography by inverting the residual anomaly after subtracting the above-mentioned contributions (Fig. 8).
The sedimentary succession was interpreted and modelled from the 3-D seismic cube generated by Alcalde et al. (2014). Eight layers were used to characterize the signal of the sediments. Three of them, corresponding to the Cenomanian-Maastrichtian, Early Cenomanian (Utrillas Fm.) and Aptian- Albian (Escucha Fm.) layers, lie above the sea level. They were used to compare the BA obtained using a reduction density of 2400 kg m −3 with that obtained using a model based on Alcalde et al. (2014). The latter map had a higher wavelength number component that needed to be filtered out, worsening the resolution of the final results and accordingly was discarded. The remaining five layers are below sea level and correspond, from shallower to deeper, to Late Jurassic to Early Cretaceous Purbeck, Dogger, marly Lias, Lias limestones and Lower Jurassic anhydrites.
The forward calculation of the sediment gravity anomaly was performed using a constant density (Table 1) for each of the five layers described above, including the topography, down to the top of the anhydrites. To solve the ambiguity generated by the unknown depth of the base of the anhydrites we used the constraints provided by the 2-D gravity models and the seismic profile shown in Fig. 5a. These data sets support the use of a constant thickness (100 m) layer of anhydrites. Furthermore, the seismic profile presented by Alcalde et al. (2014) shows a general homogeneity in the thickness of the layers for the study area. A bottom flat imaginary boundary for the Keuper rocks was used in order to avoid cutting the overlaying strata so no errors were carried into the inversion process. Finally, the long wavelength/deeper contribu- tion to the BA was removed by using the simplest first-order polynomial.
Three cross sections (Fig. 6) were extracted from the results of the gravity inversion. These include the sedimentary layers down to the anhydrites and the basement. Two of them coincide with the location of the 2-D gravimetric models, NE-SW and WNW-ESE (Figs. 3, 5a, c) and the third one has a NW-SE direction (Fig. 3). The cross sections have been used to make a comparison with the previous 2-D models shown in Fig. 5 and to assess the basement structure crossing the main tectonic features that are also highlighted by the derivatives (Fig. 4c, d).

Results from the microgravity data
The BA map (Fig. 3) shows negative values ranging between −53.2 and −48.5 mGal, which are characteristic of the Iberian Peninsula continental crust. The BA map is portrayed by a minimum located in the west-central part of the survey area where two parallel rows of low values are identified with a NW-SE direction. From there, the anomaly increases towards the south and the north, suggesting a concave ENE-WSW direction structure, with maximum values found in the SE area. The gravity gradients found to the south of the minimum are strong.
The separation of the long and short wavelength components by upward continuation of the gravity field shows that the regional anomaly (Fig. 4a) is characterized by an E-W directed low in the central/north-western area. This low shows values ranging between −52.6 and −49 mGal and is bounded by highs to the SE and partly to the NW. The residual anomaly resulting from subtracting the regional anomaly from the Bouguer anomaly map (Fig. 4b) shows variations of 1.5 mGal. The residual anomaly is characterized by NW-SE lows located in the centre of the study area and bounded by highs to the S and NW. Furthermore, an E-W strong gradient in the southern part of the map appears.
The derivatives presented in Fig. 4c and 4d aid in the interpretation of the structures affecting the Hontomín area. Above all, the gravity gradient existing to the south of the gravity data set and already evidenced in Figs. 3 and 4b is clearly imaged by both derivatives. The vertical derivative shows, as does the residual map in Fig. 4b, NW-SE minima that appear to be cross-cut by an E-W maximum in the southern part of the map. The short wavelength that these features present in the vertical derivative map indicates that they respond to shallow density contrasts.
The computed 2-D models have helped to capture the deep structure of the area. These showed to be more sensitive to the variations in the shallower layers and in the basement topography than to those in the mid-sections of the sedimentary succession. This can be due to the small differences in density found across the Jurassic succession. Model 1 has a NE-SW direction (Fig. 5a) and crosses the major faults of the area, i.e. south and east faults (Alcalde et al., 2014). Both faults are identified in the model to have affected the base-  Fig. 5a and c, (b) regional component of the BA in Fig. 3, (c) gravity residual anomaly resulting from subtracting the map in (a) to the BA map in Fig. 3 and (d) final anomaly grid resulting from subtracting the map in (b) to that in (c). This grid will be the one used for inversion procedure. ment and the sedimentary succession. While the south fault affects all the sedimentary packages, the east fault seems to be fossilized by the late Cretaceous sediments. The structural configuration of the profile shows a dome-like structure with a thick layer (up to 2000 m) of Keuper evaporites. The Jurassic succession is affected by a minor fault that affects the lower Cretaceous as well.
Model 2 has an approximate E-W direction and crosses the east fault (Fig. 5c). The vertical offset of the east fault is around 400 m and has a 45 • dip. The basement in the western sector shows a progressive uplift towards the E. The dome-like structure is recognized in this model as well as in Model 1 and is characterized by a growth of the Keuper evaporites thickness, forming a gentle fold on the overlaying Jurassic layers.
The processing workflow applied to the BA data to perform the inversion of basement topography has led to a series of maps presented in Fig. 7. Figure 7a represents the gravity response of the sedimentary succession interpreted by Alcalde et al. (2014). The forward calculation includes the Cenozoic down to an internal boundary within the Keuper rocks. Gravity minima are found to the SE of the study area and describe an irregular shape in the SE corner. These minima are in accordance with the location of a deep basement as deduced by the gravimetric models, allowing the existence a thicker Keuper layer. The maximum is located to the W and has a round shape. Another relative maximum can be found in the centre of the study area, coinciding with the deepest point of the basement in the 2-D gravimetric models. Figure 7b presents the calculated long-wavelength filter, representing the regional, and shows a gentle negative gradient towards the NW. The values have a range of less than 3 mGal. Given the reduced dimensions of the survey area, it is difficult to discern the depth of the source we are filtering out. In any case, we believe that the signal is related to density variations within the upper crust and below the top of the basement, and may be either compositional (a lithological gradient) or structural (a structure imposing a gravity gradient, e.g. low offset fault).
The resulting final grid used to perform the inversion of the gravity data is that obtained from subtracting the regional contribution in Fig. 7b to the BA map without the sediment gravity signature (Fig. 7c). This approach represents the most accurate model that we could build to ensure the best outcome of the gravity inversion, and it is shown in Fig. 7d.
The resulting final grid again shows an E-W minimum in the central part, reaching values of −2.1 mGal. Maxima appear once more in the SE corner of the area and reach values above 1.4 mGal. A strong E-W gradient exists again in the southern part of the study area. This map allows us to undertake the inversion of the geometry of the boundary, given that anomaly, i.e. the top of the basement/base of the Keuper layer.
The inversion has been performed within the Geosoft platform, using a GM-SYS 3-D module. The goal of this procedure is to obtain the geometry of the basement since the densities of the intervening lithologies are well constrained by well data. The initial set-up consisted of a two-layer model, the upper one formed by Keuper salts and the bottom one by the basement. A mean density of 2520 kg m −3 , which includes all the layers above the Keuper, even the ones located above sea level, was used as background density. The results of the inversion are shown in Fig. 8a and picture a deep basement ranging roughly from −2000 to −2650 m below sea level, with a mean depth of −2300 m. Its geometry shows a pair of minima in the central-west area, bounded at the NW and SE by local maxima. The minima lay parallel to a welldefined gradient zone that has been identified in most of the gravity maps generated for this paper. This feature seems to represent one of the main structural elements of the area and, according to its position, is interpreted as the south fault (Alcalde et al., 2014, and Fig. 5a).
The geometry of the top of the basement allows us to calculate the thickness of the Triassic Keuper salts (Fig. 8b) by calculating the space between the top of the basement and the bottom of the anhydrites layer. Results show that the Triassic Keuper appears to have a mean thickness of 1660 m with a maximum of 2020 m around the centre-east of the area. This is in agreement with the data derived from exploration well Rojas NE-1, near Hontomín, where a total thickness of 1400 m was drilled (Carola et al., 2015). The borders of the study area appear to be the location where the Triassic Keuper is thinnest, with values of about 1200 m. In general, the thickness of the layer is conditioned by the mobility of the salt and the migration pathways generated by the tectonic events.
Three cross sections (location in Fig. 3) presented in Fig. 6 were extracted from the results of the gravity inversion, which at the same time includes constraints from 3-D seismic. Two of the cross sections coincide with the two 2-D gravity models and show a relatively good correlation with them ( Fig. 5a and c). In particular, cross section 1 running in a NE-SW direction has a good correlation with the 2-D Model 1 and shows a deep basement affected by the south and east faults as well as a thick Keuper layer overlaying the basement. The third profile runs in a NW-SE direction and accounts for the NW basement high presented in Fig. 8a.
Finally, a THD has been applied to the residual gravity anomaly obtained from the Hontomín area (Fig. 7d) and is presented in Fig. 9. Results of this process are relevant for identifying gradient zones and accordingly fault zones, and will be discussed in the next section.

Discussion
The outcome of the analysis of the high-resolution BA map of the Hontomín area and its integration with the 3-D seismic results acquired over the same area (Alcalde et al., 2014) results in a coherent model of the structure of the dome that sheds some light about its evolution. Even though a good agreement between the gravity and seismic data exists regarding the topography of the basement, some issues remain unsolved in the NW sector and will be discussed further in Sect. 5.2.

Structural setting of the Hontomín dome basement
The Hontomín dome has been described as an extensional forced fold (Tavani et al., 2013) resulting from the extension and reactivation of Permo-Triassic normal faults affecting the basement. The presence of a detaching level represented by the Keuper Triassic salts produced the forced folding in the upper sedimentary cover by migration of these salts (Tavani et al., 2013;Carola et al., 2015). For the area of study, Carola et al. (2015) have suggested a thin-skinned configuration based on surface geology, well data and vintage seismic lines.
Here, we present a comprehensive model, from surface to basement, centred in the surroundings of the spot where the CO 2 storage site is currently being developed. The new model suggests a local thick-skinned deformation style for the area with two major faults affecting the basement, namely the south and the east faults, already described by Alcalde et al. (2014). This model defines an area divided into three main blocks, south, centre-north-west and north-east.
The THD presented in Fig. 9 delineates high-gradient areas affecting the residual gravity data set of Fig. 7d. These zones are interpreted as faults (discontinuous lines in Fig. 9) and are compared with the ones described by Alcalde et al. (2014). Among these, one striking ∼ E-W stands out in the southern area and is interpreted as a major fault, coinciding with the south fault as defined by Alcalde et al. (2014). The gravity gradient of this fault is observed in most of the gravity maps produced in this paper, including the residual gravity map and derivatives shown in Fig. 4 as well as those in shown in Fig. 9, indicating that the south fault affects sediments and basement. Furthermore, the cross sections derived from the 3-D model (Fig. 6) show an offset in the basement of 150 m for the south fault. Another fault parallel to the south fault is shown in Fig. 5a, close to the thickest interval of Keuper evaporites, and features a normal offset of 200 m over the basement. These two faults combined create a downward displacement of the basement of about 350 m. The south fault is thought to be a branch of the rightlateral Ubierna Fault affecting the basement and all the stratigraphic succession and conditioning the structural setting of the area. This fault has been affected by two deformation stages: a Late Jurassic-Early Cretaceous extensional stage (Tavani and Muñoz, 2012) and a later Cenozoic compression related with the formation of the Pyrenees. The kinematics and sedimentary history of this fault were interpreted by Alcalde et al. (2014). They suggest a flower-like structure associated with the strike-slip movement of the Ubierna faults. The Jurassic-Lower Cretaceous succession shows a thicker sedimentary record to the NW of the south fault, suggesting a normal displacement of the hanging wall during this period.
Another fault interpreted from the map in Fig. 9 is located to the NE, strikes NW-SE and is correlated with the east fault of Alcalde et al. (2014). This fault affects the basement and the sedimentary succession up to the Cenozoic packages, which are not affected by it. The fault has two minor associated faults (Fig. 9) that strike in the same direction but have less extension. The vertical motion of this fault was also described by Alcalde et al. (2014). Here, a downward displacement of the SW block during the Jurassic period has been assigned by the information extracted from the exploration wells. The 3-D model created after the inversion of the gravity data shows a SW dipping fault with a normal sense of motion and offsets of around 400 m, which are comparable to the offset accumulated by the southern fault. The east fault, as interpreted in the seismic data set, is not as clear in the derivatives of the BA data ( Fig. 4c and d), although it is visible and clear in the residual maps (Fig. 4b), in the basement depth maps (Fig. 8a) and in the THD performed over the final anomaly grid used for the inversion (Fig. 9). This might be due to the fact that it does not affect the entire sedimentary sequence. In general, the south and east faults seem to have been reactivated during the Jurassic-Lower Cretaceous extension, generating a sunken basement block. Figure 10. Faults interpreted from the THD of the residual gravity map and shown in Fig. 9, superimposed on the basement topography map. Figure 9 also shows good agreement existing between the south and east faults as deduced from gravity and 3-D seismic data. Good correlation exists too between minor fractures associated with these faults, like an E-W fault located to the north of the south fault and interpreted from both data sets.
Finally, the two 2-D transects modelled from the Bouguer Anomaly map (Fig. 5a and c), both of which run over one or two of the Hontomín 1 and 2 boreholes, show a similar set of faults and the same reservoir-dome-like structure. However, these models cannot be compared to the seismic lines at depth since the seismic data do not reach the basement level.
The faults interpreted in the central sunken block by Alcalde et al. (2014) are not clearly recognizable in the THD map ( Fig. 9), which shows an irregular anomaly pattern in the central domain. This could be due to either the confinement of faulting within the cover succession or to the minor offsets associated with these faults not generating strong enough gravity gradients to be recognized with this method.

Structural setting of the sedimentary cover
Another set of faults can be interpreted from the THD of the reduced BA map (Fig. 9). These strike NNW-SSE and have limited length. Even though they have a weak signature, they are similar to those interpreted from the seismic data striking in an almost N-S direction (Alcalde et al., 2014), which were described as active during the Liassic period. However, there is a discrepancy regarding their area of influence, as the latter were interpreted to be associated with a possible Triassic extensional normal faulting affecting the basement, which triggered the movement of the Triassic evaporates. Furthermore, either they do not affect the basement topography (Figs. 8a, 9 and 10), except locally to the SW, or at least they do not produce a basement offset identifiable in the gravity data. This suggests that this set of faults may only affect the Jurassic and Triassic succession and that they can be associated with the movement of the Triassic evaporites towards the basement wall, causing the dome growth.
Despite the general good agreement between the basement model presented here and that presented by Alcalde et al. (2014), there are still some differences, especially concerning the NW sector. Here, the basement topography appears to be higher when deduced from gravity data than in the seismic model. Three hypotheses can be proposed to explain this discrepancy (Fig. 11). The first (Fig. 11a) is that the basement high deduced by the gravity study is real. This basement high could be explained as a rollover anticline generated by the listric geometry of the south fault. The remain- Figure 12. Correlation of the joint gravity and seismic model with that presented by Carola et al. (2015). Note the thickening of the salt layer to the NW. A new fault is proposed to accommodate the increased thickness. The location of the fault is speculative. The tilt of the hanging-wall is taken from the south fault.
ing hypotheses imply that the gravity high observed in the NW sector is not only produced by a basement high, but also by a complex organization of the stratigraphic succession. The second hypothesis (Fig. 11b) has a reduced Keuper thickness in this area, which would imply a general increase of the gravity signal. This could be explained as the result of the migration of the Keuper evaporites towards the dome and would imply a thickening of the Jurassic-Recent stratigraphic succession in this zone. Figure 11c shows a third possible scenario where the positive gravity anomaly is related to the occurrence of a dense stratigraphic unit, which could correspond to Triassic anhydrites laterally passing to the less dense Keuper evaporites of the dome domain. Nevertheless, we cannot rule out the possibility that this structure is related to the edge effects generated after filtering a regional component in such a small area.
The basement high (2100 m b.s.l.) at the NW end of the gravimetric model also contrasts with the position of the basement top at −2800 m b.s.l. 5 km to the north as interpreted by Carola et al. (2015) based on seismic data (Fig. 12). These authors propose a thin-skinned model detached at the top of the basement and associated with a large offset on top of a low angle footwall ramp. However, the staircase geometry of the basement top across the Hontomín structure ( Fig. 12) is in our opinion in disagreement with this model, since the cover does not show a structural culmination associated with the NW basement high, which should be expected in a thin-skinned model. Besides, the monoclinal attitude of the Upper Cretaceous limestones on the southern border of the Hontomín structure above the south fault and the homogeneous uplift of this layer north of this point sug-gests that this feature is mainly associated with the inversion of the south fault and hence is a thick-skinned model. In this scenario, the basement step NW of the gravimetric model would be associated with a new normal fault (Fig. 12), which was not inverted during the compressive stage. The NW fault would be mainly responsible for the extensional forced folding occurring during the rifting stage as already proposed by Tavani et al. (2013). At a regional scale, the thickness of the Triassic salt layer in the Plataforma Burgalesa is very inhomogeneous from north to south. Exploration wells have drilled different thicknesses, ranging from 2000 to just 300 m (Carola et al., 2015). This is also the case in the Hontomín area, where the average thickness of the Keuper evaporites is 1660 m with a maximum of 2020 m and a minimum of 1200 m. This thickening generates a gravity low observed in the BA map, which appears somehow displaced from the south fault, as is also clearly pictured in the 2-D models and the cross sections derived from 3-D models (Figs. 5a, c, 6). These offsets of the gravity minima provide evidence of the mobilization of the evaporites towards the walls of the principal faults (south fault and the east fault) related to the migration pathways forced by the Cenozoic compressional stage (Tavani et al., 2013) but also to some extent to salt mobilization during the rifting stage as a consequence of the sedimentary load of the syn-rift succession on the hanging wall of the basement faults.
We have built a consistent structural model for the Hontomín CO 2 storage site from the surface down to the basement. The model relies on the integration of high-resolution 3-D gravity and seismic data. The microgravity analysis of the area has allowed us to further constrain its tectonic setting, basement geometry, faults relationships and Triassic salt thickness. The new model reveals a thick-skinned tectonic setting configured by three sets of faults, all of them correlatable with those identified by Alcalde et al. (2014). Two of them, namely the south fault and the east fault, striking ENE-WSW and NW-SE respectively, clearly affect the basement. The south fault also affects the entire stratigraphic succession above it, while the east fault affects the basement and the sedimentary succession up to the Cenozoic sediments. We propose a halfgraben-like structure for the configuration of the south fault. Accordingly, a basement high located in the NW section of the area and not affected by any structure might suggest that the south fault has acted as a listric fault, configuring the basement topography in the NW sector. A third set of faults striking NNW-SSE has been identified. However, its gravity signature in the residual gravity map is not very conspicuous, suggesting that they may affect just the Triassic Keuper evaporites and the Jurassic succession above but not the basement. Accordingly, they have been interpreted to be the result of the movement of the salt towards the main faults. However, it is also possible that they affected the basement but produced small offsets, as previously suggested by the 3-D seismic data. This paper demonstrates that the integration of high-resolution geophysical data sets is a magnificent tool for unravelling the structure of geologically complex areas.