Edinburgh Research Explorer 3-D reflection seismic imaging of the Hontomín structure in the Basque–Cantabrian Basin (Spain)

. The Basque–Cantabrian Basin of the northern Iberia Peninsula constitutes a unique example of a major deformation system, featuring a dome structure developed by extensional tectonics followed by compressional reactivation. The occurrence of natural resources in the area and the possibility of establishing a geological storage site for carbon dioxide motivated the acquisition of a 3-D seismic reﬂection survey in 2010, centered on the Jurassic Hontomín dome. The objectives of this survey were to obtain a geological model of the overall structure and to establish a baseline model for a possible geological CO 2 storage site. The 36 km 2 survey included approximately 5000 mixed (Vibro-seis and explosives) source points recorded with a 25 m inline source and receiver spacing. The target reservoir is a saline aquifer, at approximately 1450 m depth, encased and sealed by carbonate formations. Acquisition and processing parameters were inﬂuenced by the rough topography and relatively complex geology. A strong near-surface velocity inversion is evident in the data, affecting the quality of the data. The resulting 3-D image provides constraints on the key features of the geologic model. The Hontomín structure is interpreted to consist of an approximately 10 7 m 2 large elongated dome with two major (W–E and NW–SE) striking faults bounding it. Preliminary capacity estimates indicate that about 1.2 Gt of CO 2 can be stored in the target reservoir.


Introduction
The Hontomín structure has been classified as an example of a forced fold-related dome structure (Tavani et al., 2013). This type of folding is generally closely related to the development of extensional fault systems (Schlische, 1995;Cosgrove and Ameen, 2000;Tavani et al., 2011Tavani et al., , 2013. They can develop under varying conditions (Braun et al., 1993;Jin et al., 2006) and their presence has been proposed in the Mesozoic Basque-Cantabrian Basin of northern Iberia (Soto et al., 2011;Tavani et al., 2011Tavani et al., , 2013. The formation prerequisite for these structures involves the existence of a ductile layer which acts as a major decoupling zone that is able to separate the deformation of the sedimentary cover from the deformation and faulting occurring at the basement. The decoupling level in this area is formed by Triassic evaporites, which separate the pre-rift Upper Triassic-Middle Jurassic folded sediments from the Paleozoic-Lower Triassic faulted rocks of the basement. The structure appears to have been generated by reactivation of the pre-existing extensional basin architecture during a compressional phase associated with Cenozoic tectonics affecting the Pyrenees. During this compressive stage the overlying Mesozoic sequence overrides the Cenozoic Ebro foreland Basin (Martínez-Torres, 1993). Thus, the area constitutes a unique case study of a major deformation system within the Basque-Cantabrian Basin in the northern Iberian Peninsula.

J. Alcalde et al.: 3-D reflection seismic imaging
The Hontomín structure, along with the Ayoluengo structure, has been targeted for hydrocarbon exploration in the northern Iberian Peninsula (Alvarez, 1994;Merten, 2006;Beroiz and Permanyer, 2011). Thus, the area has been extensively studied since the early 1970s. A number of 2-D seismic reflection transects were acquired with the aim to unravel the regional structure (Beroiz and Permanyer, 2011). The subsurface stratigraphy includes a carbonate Jurassic sequence that hosts a deep saline aquifer and a high-quality seal formation. This stratigraphy motivated new interest in the area as the structure could be considered feasible for geological storage of CO 2 (Orr, 2004;Baines and Worden, 2004;Alcalde et al., 2013).
Geophysical methods have been broadly used for the characterization of geological structures (Martí et al., 2002;Preston et al., 2005;Förster et al., 2006). Amongst these, seismic methods provide high-resolution images and baseline models of hydrocarbon producing and/or storage reservoirs (Juhlin et al., 2007;White, 2009). The interest in the Hontomín structure motivated the acquisition of a number of multidisciplinary data sets intended to characterize the Hontomín site (e.g., Alcalde et al., Elío et al., 2013Ogaya et al., 2013;Ugalde et al., 2013). The active seismic experiments carried out in the area were initially discussed in Alcalde et al. (2013). These included a 3-D seismic reflection data set, acquired with the aim of obtaining a 3-D image of the Hontomín structure. Constraining the tectonic framework of the area and providing a geological and seismic baseline model were the main aims of the survey.
The study area is geologically complex and, at the surface, topographically irregular. One of the difficulties encountered during processing is the existence of velocity inversions at relatively shallow depths which complicate the determination of optimum static corrections. In general terms, the Hontomín area comprises a very heterogeneous folded structure, composed of mixed sedimentary marine carbonate and continental siliciclastic sediments. According to well-log data, the reservoir-seal system within the Hontomín structure is a carbonate sequence located at approximately 1500 m depth. Reflection seismic methods often encounter difficulties when attempting to image carbonate reservoirs. Suggested reasons behind this are the absence of distinct layering during the deposition and diagenetic alterations of the carbonate sediments (Masaferro et al., 2003;Braaksma et al., 2006;Von Hartmann et al., 2012), abrupt impedance and velocity changes (Rudolph et al., 1989;Janson et al., 2007), and more generally, high structural complexity (Phipps 1989;Rudolph et al., 1989;Erlich et al., 1990;Kenter et al., 2001). This was an important issue during velocity analysis. Other problems encountered related to acquisition and processing are presented and discussed in this work.
The revised processing scheme used in the present study was adapted to the characteristics of the Hontomín subsurface. The resulting final image features higher coherency than that presented in Alcalde et al. (2013) throughout the entire seismic volume. The improvement is mainly due to the calculation of a more precise static model, the improved noise reduction and the source wavelet matching method applied. This enables a better analysis and identification of the subsurface structures, providing a definitive seismic reference model of the study area. Although the focus of this paper is the acquisition and processing of the 3-D data set, a preliminary interpretation is also presented.

Geological setting
The Hontomín structure ( Fig. 1) is bordered by the NW-SE Ubierna fault to the southwest (Tavani et al., 2011), the Poza de la Sal diapir to the north (Quintà et al., 2012) and the Ebro and Duero basins to the east and southwest, respectively. The study area is located within in the North Castilian Platform, a Mesozoic sub-unit on the western part of the Basque-Cantabrian basin (Serrano and Martinez del Olmo, 1990) that developed during the opening of the Atlantic Ocean and the Bay of Biscay (Ziegler, 1989). This Mesozoic extensional regime generated a synformal domain in the region (Pujalte et al., 2004;Quintà and Tavani, 2012), characterized by a thick sedimentary sequence (Quesada et al., 2005). Due to the subsequent Alpine compressional regime, small-scale inversion structures were produced, leading to reverse, rightlateral and left-lateral reactivation of the pre-existing faults (Tavani et al., 2011;Quintà and Tavani, 2012;Tavani et al., 2013).
The extensive exploration for natural resources within the area (Alvarez, 1994;Merten, 2006;Beroiz and Permanyer, 2011) produced a relatively large amount of geophysical and geological data. This includes vintage 2-D seismic reflection data, well-logs, borehole core samples, etc. The existing data indicate that the reservoir and seal formations are Jurassic in age, and form a slightly elongated dome-like structure with an overall aerial extent of 5 km × 3 km. The area features low natural seismicity as stated in Ugalde et al., (2013). From a stratigraphic perspective, the lowest Mesozoic succession (Fig. 1b) is the Keuper facies, which forms the core of the target dome. The "Carniolas" unit is composed of evaporites, dolomites and marls (sabkha facies) and lies over the Keuper facies, representing the boundary between the Triassic and Lower Jurassic (Quesada et al., 2005;Pujalte et al., 2004). The Lower Jurassic sediments can be described with four units: Limestone Lias and Marly Lias, and Pelletic dogger and Limestone dogger. They are characterized by shallow marine carbonate ramp sediments and hemipelagic ramp sediments. The Upper Jurassic-Lower Cretaceous transition is constituted by the Purbeck facies, formed by clays, sandstone and carbonate rocks. The remainder of the Lower Cretaceous succession consists of siliciclastic sediments from Weald, Escucha and Utrillas facies, mainly composed of a series of sandstones and shales, and occasionally microconglomerates. The uppermost rocks (outcropping in the  Hontomín area) are Upper Cretaceous carbonates and Cenozoic rocks (lacustrine and detritic) lying unconformably over the Mesozoic successions. The structure is bounded to the south by a strike-slip fault associated to the Ubierna Fault (Tavani et al., 2011).
Four pre-existing boreholes (H1, H2, H3 and H4 in Fig. 2), reaching depths between 900 and 1832 m, sampled the Mesozoic succession. These boreholes were drilled between the 1970s and the 1990s for oil exploration purposes and some were still producing in 2009. Beneath the producing level, a saline aquifer formation is present. It is located at approximately 1500 m depth, within Early Jurassic (Raethian-Sinemurian) shallow marine carbonate ramp sediments. This potential CO 2 storage formation is composed of a dolostone unit (upper part of the Carniolas unit), overlain by an oolitic limestone formation (Limestone Lias), with an average thickness of approximately 80 m. This reservoir level is fairly porous (between 9 and 17 %, Ogaya et al., 2013) and contains saline water with more than 20 g l −1 of NaCl. The porosity of the Carniolas is believed to be the result of secondary dolomitization and different fracturing events. There is no available information on its permeability. The potential upper seal unit is comprised of marlstones and black shales from the hemipelagic ramp, corresponding to the Marly Lias unit (Pliensbachian and Toarcian in age).

Seismic data acquisition
The acquisition parameters (Table 1) were selected based on the knowledge derived from previous studies in the area (Gemodels/UB 2011, internal report). These previous studies included surface geological mapping, reprocessing of vintage 2-D seismic exploration surveys and correlation with the  exploration boreholes, which provided overall information of the subsurface structures. The survey size was designed to cover the entire Hontomín structure.
The 3-D seismic survey geometry included 22 source lines perpendicular to 22 receiver lines, with intervals of 25 m between sources and receivers (Fig. 2). The distance between receiver lines was 250 m, and source lines 275 m. The survey was acquired in five stages corresponding to different patches (swaths), with patches at the edge of the survey consisting of six active receiver lines, and the inner patches consisting of 10 receiver lines. A maximum of 120 channels were active per receiver line, resulting in a maximum of 1200 traces per shot gather.
The survey was originally designed to be fully acquired with a Vibroseis source (M22 vibrators). However, logistical issues (mainly due to the rough topography in the eastern part) forced the use of explosives in about 24 % of the acquired data. The selected parameters for the Vibroseis source were two 16 s sweeps, with an 8-80 Hz bandwidth. The explosive source consisted of 450 g of explosives, deployed in three 1.5 m deep holes for each shot point. The total number of source points was 4818, resulting in 4 333 605 traces.
The use of two different sources during the acquisition resulted in different seismic signatures in the shot gathers depending upon what type of source was used (Fig. 3). The explosive shot records are more strongly contaminated with ground roll and low-frequency guided waves than the Vibroseis shot records. The frequency content of the Vibro-seis shots was determined by the sweep bandwidth of the source (i.e., 8-80 Hz), whereas the explosive data, which feature lower amplitude, resulted in a broader bandwidth, up to 100 Hz (Fig. 4).

Seismic data processing
The relatively regular pattern of the acquisition geometry (Fig. 2) gives a natural bin size of 12.5 m × 12.5 m for the CDP bins, resulting in a maximum CDP fold of 36 traces/CDP (Fig. 5a). This fold value is relatively conventional in land surveys and for reservoir characterization (Ashton et al., 1994;Roche, 1997) and affordable from the academic point of view. According to previous information of the study area (Geomodels/UB 2011, internal report; Ogaya et al., 2013), dipping events steeper than 16 • are rather unexpected, and therefore, quasi-horizontal dips were assumed in the processing. Source-receiver offsets of up to 2500 m ( Fig. 5b) were enough to constrain reflections from a depth of 1500 m (1150 ms two-way travel time (twt), assuming an average velocity of 3500 m s −1 ). Furthermore, at travel times corresponding to this target depth, surface-wave energy is considerably decreased. The swath acquisition provided a well-sampled azimuthal coverage (Fig. 5c), which assures the three-dimensionality of the data. The available well-log data indicate that the lithologies with higher impedance contrast are located above 1400 ms twt in most cases. Taking this into account the data were processed down to 1500 ms. The processing flow (Table 2) was designed and implemented in order to enhance the S/N ratio and image the target structures. Alcalde et al. (2013) applied a standard processing to the 3-D data set. A review of this preliminary processing was conducted to identify some key data/noise issues associated with the data set. These issues included the identification of the noise sources present during the acquisition, the shadow zone related to a near-surface velocity inversion and the source wavelet matching issue. The new processing workflow presented in this paper and described below addresses these issues and takes a more detailed look at the signal processing in order to obtain a final image suitable for subsequent detailed interpretation.

Source matching
The use of multiple sources in seismic data acquisition introduces additional complexity to the processing, due to the need to homogenize their phase characteristics (Steer et al 1996;Feroci et al 2000). Since both Vibroseis and dynamite sources were utilized in the survey (Fig. 2), efforts were made to match source wavelets and bulk time shifts between the dynamite and Vibroseis traces, in order to improve the lateral continuity of the final stacked volume. Although the conventional approach is to convert the Vibroseis source data to minimum phase, given that most of the survey consisted of  Vibroseis shots, with the dynamite shots confined to the survey edges, the minimum phase dynamite wavelet was converted to zero-phase to match the cross-correlated Vibroseis wavelet. In order to do this, two wavelet matching processes were tested: (1) applying trace wide time and phase shifts and (2) a convolutional match filter and time shift. The match filter used to convert the dynamite wavelets to zero-phase was based on a wavelet extracted from 10 dynamite shots from different parts of the survey.
In order to assess the best parameters for each different wavelet matching method, co-located dynamite and Vibroseis trace pairs were selected. For the time and phase-shift wavelet matching method, a range of different values of both were applied to the dynamite trace. Thus, the correlation coefficient between the two traces in a given time window was assessed. This analysis was also performed in the same time window after applying the match filter and variable time shift. The optimum parameters were selected based on those which gave the highest correlation coefficient. It was observed that, based on tests performed with synthetic data, the effects of environmental noise could be overcome by stacking the correlation coefficients across many dynamite-Vibroseis trace pairs. Figure 6 shows the stacked correlation coefficients for 1080 pairs of dynamite and Vibroseis traces, located within 25 m of each other, taken from 10 different shots. The best match between the dynamite and the Vibroseis traces is obtained at time and phase shifts of 12 ms and −150 • ,

Fig. 4.
Frequency spectra of the shot gathers shown in Fig. 3: shot gather 1820, acquired with Vibroseis source (dark grey) and shot gather 1158, acquired with explosive source (light grey). The highcut filter in the spectrum of the Vibroseis data was applied in the field.
respectively. For the match filter method, the optimum time shift to use in addition to the match filter was 24 ms. Certain corrections were performed in order to diminish the disparity in the amplitudes and frequency content of the shot-gathers (Fig. 4). The frequency contents were equalized  by the use of combined bandpass filtering and spectral equalization in the pre-and post-stack stages. The high-amplitude contrast forced the use of AGC and trace balancing (250 and 300-1300 ms time windows, respectively), since the spherical divergence correction did not completely equalize the amplitudes.

Static corrections
Elevation and refraction statics are key processes for obtaining a high-quality final seismic image (Juhlin et al., 2007;Malehmir and Bellefleur, 2009;Malehmir and Juhlin, 2010). In our case, they were essential for improving the final image due to the large topographic changes (up to 200 m) and heterogeneous geology present in the study area. The entire data set was referenced to a datum of 1070 m a.m.s.l., just above the highest altitude of the survey. The replacement velocity, based upon apparent velocities in the shotrecords, was 3500 m s −1 . The resultant elevation static model adjusted the data by 67 ms on average, with maximum corrections of 10 ms.
The refraction static modeling was an arduous process due to the heterogeneous characteristics of the subsurface. One of the most interesting features observed is the lack of clear first arrivals in the data at offsets greater than 300 m (black arrows in Fig. 3c, d). This signal disappearance is observed as a sharp (black arrows in Fig. 3c) or as gradual loss of amplitude (black arrows in Fig. 3d) with increasing offset. This phenomenon is observed in the whole survey for both source types. At far offsets, this results in a lack of signal in the upper 500 ms (Fig. 3e), i.e., a "shadow zone". The existence of this shadow zone added difficulty to the pick-ing of first arrivals in the near-offset traces, and impossible at offsets larger than about 700 m. In spite of this difficulty, more than 670 000 first breaks, ranging from 0 to 500 m offsets, were manually picked. The refraction statics were calculated using a two-layer model, with an input-fixed upper-layer velocity (at 1070 m m.s.l.) of 1900 m s −1 , and a lower half-space velocity (at a starting depth of 40 m below the surface) of 3500 m s −1 . The resulting model refraction statics ranged from 0 to 40 ms (RMS of approximately 8 ms), whereas the residual refraction statics were always below 4 ms. Static corrections, particularly the refraction static corrections, increased the quality of the final image and improved the coherency of the reflections in the data (Figs. 7c and 8c).

Stack, post-stack processing and migration
Velocity analysis after static corrections provided further support for the existence of a velocity inversion between 200 to 500 ms. The NMO-velocity reduction in this interval ranged from 300 to 600 ms. The CDP gathers show disrupted reflections before and after NMO corrections (Fig. 9). An example of a stacked section (up to step 14 in Table 2) is shown in Fig. 10a. After a satisfactory NMO model had been obtained, a 3-D Kirchhoff DMO code (based on Hale and Artley, 1993) was applied to the pre-stack data, using the NMO velocities as initial input (Fig. 10b). Although the Hontomín structures are not steeply dipping, the 3-D DMO processed volume shows increased coherency and improved definition of the flanks of the dome, as well as improvement in the continuity of reflections, as compared to the stacked data without DMO applied (black arrows in Fig. 10a, b). Post-stack F XYdeconvolution in the inline and crossline directions was implemented to remove random noise and further increase coherency. Different finite-difference and phase-shift time migration codes were tested. We chose a 45 • finite-difference post-stack algorithm, due to the spatial variation of the velocities present in the study area.

Analysis of the seismic image
A section comparing the preliminary image obtained by Alcalde et al. (2013), and the same section obtained in this work can be found in Fig. 11. In the revised work flow, the subsurface is well imaged in the migrated volume from approximately 75 ms down to 1350 ms, completely covering the target structure. Significant differences between the two volumes can be observed. In general, the revised work flow exhibits less low-frequency background noise and an overall increase of the signal-to-noise ratio. The northern part shows clearer reflections and improved identifiability of the defining structures. The revised work flow also has lower  Fig. 4) comparing the amplitudes and frequencies of the shot 1820 (Vibroseis source, dark grey) and 1158 (explosive source, light grey). Note that frequency graphs are shown only for qualitative comparison (values of the frequency graphs are the same as in Fig. 4, not shown here for display purposes). 711 Fig. 7. Closeup of the shot gather 1820 (see also Fig. 3) after different steps in the data processing: (a) raw shot; (b) same as (a) plus spherical divergence and surface consistent deconvolution; (c) same as (b) plus elevation and refraction statics; and (d) same as (c) plus time-variant frequency filter and spectral whitening. Each plot also contains a frequency spectrum graph (similar to that of Fig. 4) comparing the amplitudes and frequencies of the shot 1820 (Vibroseis source, dark grey) and 1158 (explosive source, light grey). Note that frequency graphs are shown only for qualitative comparison (values of the frequency graphs are the same as in Fig. 4, not shown here for display purposes).

FIG. 9:
Example of two CDP-gathers from different areas before and after NMOcorrection, and the two NMO velocity models applied in each case. The highlighted boxes mark a reflectivity area whose reflections are not continuous even after a correct NMO correction. Fig. 9. Example of two CDP gathers from different areas before and after NMO correction, and the two NMO velocity models applied in each case. The highlighted boxes mark a reflectivity area whose reflections are not continuous even after a correct NMO correction. migration-related artifacts, especially in the edges of the section.
The lateral continuity of the reflections is limited in the whole seismic volume (e.g., sets E and F in Fig. 12). This could be due to the influence of the shallow velocity inversion, as well as to the existence of heterogeneities associated with small-scale fracture zones. However, the quality of the processed data is enough to allow us to perform a detailed analysis of key seismic events/facies and their distribution within the migrated seismic volume (Figs. 12a and Fig. 13). Line α-α (location in Fig. 2) crossing boreholes H4, H2 and H1 is used in this section to analyze the seismic image and correlate it with the main geological units (Fig. 12). This characterization is applicable to the whole seismic volume.
Three domains were identified by the analysis of the seismic facies in the seismic data set. The first domain generally contains low-amplitude, discontinuous reflections (A). It is located at the bottom of the seismic volume, below approximately 1100 ms. It shows a reduction in amplitudes compared to the interval immediately above, and less continuity in the reflections. Within this domain, a brighter set of parallel reflections (B) is found between inlines 1305 and 1340. Overlying the A and B facies sets, the second domain is composed of a sequence of more continuous, high-amplitude reflections (C). They show an asymmetric dome-shaped geometry with an average thickness of 300 ms (Fig. 12). The boundary between the A and C reflection packages is chaotic across most of the seismic volume. The C set is characterized by higher amplitudes compared to the surrounding seismic facies, with dominant frequencies around 30 to 40 Hz. Its amplitudes decrease towards the end of the second domain.
Over the dome flanks, the third domain begins with a set of reflections (D) overlying the C set. Set D is characterized by very high amplitudes and a 80-120 ms −1 thickness. However, the D package is only observed in the northern portion of the study area. The overlying E set is characterized by variable amplitude and non-planar reflections. A high-amplitude reflection at approximately 575 ms bounds the E set. Above this, a set of higher amplitudes and more continuous parallel reflections appear (F). This package has an approximately constant thickness of 400 ms. In Fig. 12, set F shows an important lateral change of the seismic properties (at inline 1260, approximately), where a more continuous character to the north becomes less continuous to the south. At the top of the seismic volume, another set of discontinuous reflections (G), nearly 200 ms thick, covers the whole area and generally extends to the surface. In the southern part, a set of southdipping reflections (H) with high-coherency overlies set G near the surface.

Correlation between seismic data and geology
The subsurface structure of the Hontomín CO 2 storage site described below is based on the seismic data, information from borehole data and previously presented information of the underground and surface geology. The available welllog data lack velocity and/or density information in approximately the first 400 m, and totally lack check-shot information. This makes it extremely difficult to obtain a precise seismic to well tie. An example of the calculation of a synthetic trace for the H3 well using the available well-log information is shown in Fig. 14 Table 2) from α-α section (location in Fig. 2 (Fig. 12b).
The Triassic package (Fig. 12b) generally has an inhomogeneous seismic pattern. It includes reflection sets A and B. The presence of a variable quantity of anhydrite layers, together with the associated increase in velocity of these layers, could explain the lack of clear seismic events. Only in spe-cific locations, such as B (Fig. 12a), can clear sub-horizontal layering be defined by high-amplitude reflections. Surface and borehole data show a non-uniform, gradual change from anhydrite to carbonate. This could be linked to the diffuse seismic response observed near the surface. The overlying Jurassic unit (Fig. 12b) is characterized by a high-amplitude reflective body observed from 750 to 1200 ms, corresponding to seismic set C. It includes the target aquifer and an overlying seal layer. The Jurassic unit is overlain by the Lower Cretaceous unit, which includes seismic sets D, E and F. Set D probably corresponds to a transitional sedimentary package, which could explain why it only appears in the northern portion of the seismic volume. The upper part of the Lower Cretaceous is formed by the Weald, Utrillas and Escucha Fm sediments (Fig. 1). Their boundaries are unclear, but they all feature high-amplitudes and gently dipping planeparallel layering. The shadow zone is probably contained within the boundary between sets F and G (150 ms TWT, approximately). This boundary divides the Lower and Upper Cretaceous. It is especially visible in the areas where the shadow zone is more noticeable (e.g., inline 1200 to 1250 in Fig. 12a). The Upper Cretaceous carbonates outcrop across the entire study area except on its southern and eastern ends, where the Upper Cretaceous beds drop and are overlain by the Cenozoic unit (set H). This package exhibits high coherency and an approximate dip of the reflections of 16 • to the south. In the eastern part, this package is subhorizontal, displaying less coherency than in the southern end.
Three main faults were identified in the seismic data set: the southern, the eastern and the central faults (Figs. 12, 13  and 15). The Southern fault is an approximately E-W striking fault zone that crosses the entire data set. It is not purely vertical, and shows some bifurcation. The Eastern fault is a major NNW-SSE reverse fault that abuts the Southern fault. A third set of small-scale faults are observed only in a few sections. One of them is the Central fault, a reverse fault mapped only in the α-α section (Fig. 12b).   Table 2. Location of the β-β section in Fig. 2.

Data acquisition and processing
Three-dimensional seismic experiments, especially in complex environments, generally require a reasonable fold in order to obtain high-quality images of the subsurface structures. The acquisition approach used in Hontomín (Fig. 2) assured that the target could be properly imaged, with long offsets and 3-D sampling across the whole survey area (Fig. 5) and with a reasonable fold (Ashton, et al., 1994;Roche, 1997). However, it was unexpected to find a sharp velocity inversion so close to the surface. The velocity inversion is associated with the Upper-Lower Cretaceous boundary (i.e., the transition from high-velocity carbonates to low-velocity, water-saturated siliciclastic sediments). According to Fermat's principle, the seismic energy reaching the velocity inversion point tends to travel preferentially through deeper (higher velocity) formations, avoiding the shallower (lower velocity) layers. This results in a complete loss of amplitudes that severely affected the quality of the data by reducing the information in the traces corresponding to the position of the shadow zone at offsets larger than 500 m (Fig. 3). Besides, the deep reflections show a severe lack of coherency across the data set (Fig. 9). The conventional NMO correction processing thus fails to some extent since reflective zones are correctly aligned horizontally, but single reflections are disrupted and the signals do not sum constructively in the stacking process (Fig. 9). This feature could be a product of the interaction of the backscattered shear-wave energy with  the deep reflections, although further studies should be performed to elucidate this aspect.
Additional difficulties arose during the data processing due to physical property changes related to the carbonate composition of the subsurface. In this sense, the carbonate sediments behave similarly to crystalline media, containing sharp velocity variations which produce a blurred seismic signature. As observed in similarly difficult reservoir imaging experiments, the most relevant processing steps for reaching a good final image were the elevation and refraction static corrections (e.g., Juhlin 1995;Malehmir et al., 2012). Although first break picking was hindered by the "shadow zone" (Fig. 3), the resultant static model was very successful in improving the seismic image (Fig. 8), which is in accordance with the heterogeneous surface geology. However, there is still room for improvement in the statics calculation. New information about the near-surface velocity field (e.g., sonic-log data, shallow refraction seismic data, surface wave tomography etc.) could help to constrain and refine the nearsurface velocity model.
For the source wavelet matching, the two wavelet matching methods (match filtering and phase and time shift) were applied to a subset of the data which was then stacked and qualitatively assessed. Both methods improved the data and when compared to each other, both were very similar in overall performance. As the phase and time shift method was simpler to apply, this method was chosen and applied to the data in the final workflow. This resulted in enhanced reflector continuity and strength in areas of the stacked section, which contained both dynamite, and Vibroseis traces in the Sonic and density logs of the H3 well in depth, with their interpreted age. A zero offset synthetic trace, calculated assuming a Ricker wavelet with a peak frequency of 30 Hz is shown (yellow). The time-depth relationship for the synthetic was set manually, due to the lack of check shot data and velocity information of the shallow subsurface. The synthetic is overlain on a portion of the β-β' section (location in Fig. 2). The Jurassic-Lower Cretaceous boundary is correlated with a dotted line. Main seismic events are also labeled in the seismic section. 720 Fig. 14. Sonic and density logs of the H3 well in depth, with their interpreted age. A zero offset synthetic trace, calculated assuming a Ricker wavelet with a peak frequency of 30 Hz is shown (yellow). The time-depth relationship for the synthetic was set manually, due to the lack of check shot data and velocity information of the shallow subsurface. The synthetic is overlain on a portion of the β-β section (location in Fig. 2). The Jurassic-Lower Cretaceous boundary is correlated with a dotted line. Main seismic events are also labeled in the seismic section.
pre-stack gathers (Fig. 8). Due to the unequal distribution of the explosive source points throughout the Hontomín survey (Fig. 2), the source matching effects are only observed at the edges of the survey, where the explosive source concentration is higher.
The use of AGC within the workflow helped to obtain a good quality final image. However, it does not preserve the original amplitude character of the data set. This could jeopardize the use of amplitude variation methods for monitoring. Modifications to the processing flow should therefore be made before this data set can be used as a baseline to a repeat survey as part of a time-lapse seismic monitoring study.

Geometry of the deep structures
In general, the work flow outlined in Table 2 provides an improved image compared to the previous preliminary image presented in Alcalde et al. (2013) (Fig. 11). Our final processed seismic volume allows an interpretation of the five main geological units present in the study area (Fig. 12). The target Jurassic dome structure has been mapped from 650 ms to 1200 ms twt (Fig. 15a). The horizons corresponding to the top and the base of this structure have been mapped throughout the volume (Fig. 13). The top surface appears first just below 1050 ms and reaches its crest above 650 ms. The base surface appears at 1350 ms in the northern part of the data set, and peaks at 880 ms TWT. They define an elongated dome like structure. The geometry of the dome is almost symmetric up to 850 ms. From this point, the dome develops into an  elongated shape in the NW-SE direction. It maintains this geometry until ∼ 725 ms. The crest of the dome (at 650 ms) is found at the eastern side of the H2 well, in an area between inlines 1270-1300 and crosslines 1130-1160 (Fig. 15). The Jurassic dome resembles a pseudo-triangle in plan view, with the two main faults, eastern and southern faults, representing two of the sides. The third side was determined subjectively at the two-way time that the dome changes from symmetric to elongated (i.e., approximately 850 ms). The calculated area of the dome is 9.4 × 10 6 m 2 .
The Hontomín structure is faulted at different scales. The southern and Eastern fault mark the southern and eastern limit of the dome, respectively. The Southern fault crosses the entire processed volume, separating the southern third of the study area from the remainder (Fig. 12 and Fig. 13). It has been related to the Ubierna fault (Tavani et al., 2013) and described as a dextral strike-slip fault that outcrops within the study area. It is easily recognizable because it coincides with a change in the dip of the reflections (i.e., increased dips to the south side of the fault coupled with a loss of lateral continuity of the events). It is difficult to assess whether a vertical offset across it is present or not, since the fault cuts through dome-shape structures and different sections give varying interpretations of the offset. The fault also contains branches that are only observable in certain sections of the seismic cube (Fig. 13). This fault could be an important conduit for fluid circulation, as observed by Ogaya et al. (2013). In addition, hydrochemistry studies in fluids derived from surface springs conducted by Nisi et al. (2013) in the eastern part of the Southern fault suggest that these fluids could have a deep origin, supporting the idea of the fluid circulation and, thus, a leaky fault. The Eastern fault presents the largest vertical offset of the study area (up to 380 m, Fig. 15b). It is nearly vertical near the surface and slightly dipping at depth. It has been mapped up to the Lower-Upper Cretaceous boundary. It marks a sharp seismic facies contrast. Other minor faults (e.g., the Central fault, Fig. 12b) with small offsets (less than ∼ 20 ms) have been identified, but not included in this preliminary interpretation.

Internal architecture
The frequency range at the target time (700-1200 ms) is from 30 to 70 Hz. This frequency range is wide enough to ensure a reasonable resolution at the aquifer level. The dominant frequency of the data ranges from 20 to 75 Hz. Thus, assuming an average velocity of 3900 m s −1 , the vertical resolution ranges between 13 to 50 m, based on the 1/4 wavelength criterion (Widess 1973;Sheriff and Geldart 1995). Based on the Fresnel zone criterion (Sheriff and Geldart 1995;Yilmaz 2001), the lateral resolution for the deepest structures at depths of 1200 and 1500 m (approximately) range from 180 to 320 m. Migration tends to collapse the Fresnel zone (and therefore increase the spatial resolution) to approximately the dominant wavelength (Yilmaz, 2001). However, the presence of noise and migration-velocity related errors could degrade the lateral resolution significantly.
A limited definition of the internal architecture was achieved in the interpretation stage. First, the carbonate units show subtle acoustic impedance changes (determined from the existing logs) and significant scattering. This makes the identification of the internal structures very difficult. Second, to date, the available velocity information from well logs is incomplete for the first 400 m. The latter issue does not allow us to tie the seismic data to the available boreholes. In this respect, new boreholes geared to hydrogeological modeling should cover this lack of information; nevertheless, these are not currently available. New processing steps should be performed in order to carry out the detailed geological interpre-tation of the 3-D volume, including a time to depth conversion of the data with the velocity information derived from logging tools, wave-equation datuming and, if possible, prestack depth migration.
In spite of this lack of internal definition, we have calculated a theoretical CO 2 capacity of the reservoir unit within the Jurassic structure. We will not take in consideration the possible interactions of the CO 2 with the reservoir formation, since this is outside of the scope of this article. Further information about geochemical and rock physics interactions, conducted inside the Hontomín project, can be found in García-Ríos et al. (2013) and Canal et al. (2012), respectively. The outcrop data indicates that Jurassic units are relatively homogeneous and their thicknesses are relatively constant. The calculation takes into account the characteristics of the top of the interpreted dome, and extrapolates it downwards. Assuming an average thickness and effective porosity of the reservoir unit of 80 m and 8.5 % (obtained from the well-log data), respectively, the calculated total pore volume in the dome reservoir is 1.6 × 10 6 m 3 . At the expected average reservoir conditions (41 • C and 15.3 MPa), the density of CO 2 is 745 558 kg m −3 . Thus, a maximum CO 2 storage capacity of 1.2 Gt of CO 2 is expected in the Hontomín Jurassic structure. This value is a rough estimate, but provides an overall idea of the scale and potential of the Hontomín site for CO 2 storage.

Conclusions
The data processing applied to the Hontomín 3-D seismic reflection data set allowed a preliminary interpretation of the main underground structure of the Hontomín area to be made. The data also characterize the Jurassic structure, which contains a potential reservoir-seal system for the storage of CO 2 in a deep saline aquifer. The structure includes a nearsurface layer characterized by a velocity inversion. This generated unexpected complications in the static correction calculations, and reduced the overall quality of the data. In spite of this, the processing steps applied to the data resulted in relatively good images that show an elongated, faulted dome structure. Static corrections proved to be the most important step within the data processing flow. Other key processing steps included the DMO correction, frequency filtering and amplitude equalization. The two source wavelets used in the acquisition were matched using a self-developed code that gave similar results to that of the standard match filtering, with an easier implementation.
Five main geological units were interpreted from the migrated volume, including the target dome corresponding to the Jurassic package. The shape of the target structure is an asymmetric dome, with a steeper flank facing southwest. Two major faults, the Southern fault and the Eastern fault have been successfully interpreted as crossing the entire volume W-E and NNW. The total area of the dome structure is 9.4 × 10 6 m 2 . The calculated pore volume for the reservoir structures within the Hontomín dome is 1.6 × 10 6 m 3 , giving a maximum theoretical storage capacity of 1.2 Gt of CO 2 . A refined interpretation of internal structures and small-scale fracture zones, as well as time to depth conversion using the forthcoming well-log information are the next logical steps to be carried out. The processed data set represents the 3-D seismic model that could be used as a baseline model within a CO 2 storage scenario.
Acknowledgements. The authors sincerely thank the Editor Hans Thybo, Peter Bergmann and the anonymous reviewer for their useful comments. Funding for this Project has been partially provided by the Spanish Ministry of Industry, Tourism and Trade, through the CIUDEN-CSIC-Inst. Jaume Almera agreement (Characterization, Development and Validation of Seismic Techniques applied to CO 2 Geological Storage Sites) and by the European Union through the Technology Demonstration Plant of Compostilla OXYCFB300 project (European Energy Programme for 534 Recovery). Additional support has been provided by Spanish Ministry of Education Science (CSD2006-00041), Generalitat de Catalunya (2009SGR006) and CSIC JAE-Doc postdoctoral research contract (E.S.). The sole responsibility of this publication lies with the authors. The European Union is not responsible for any use that may be made of the information contained herein. Juan Alcalde is being currently supported by the Fundación Ciudad de la Energía (CIUDEN) research training program. We sincerely thank CGGVeritas for their assistance during the acquisition. Processing has been carried out using GLOBE ClaritasTM, under license from the Institute of Geological and Nuclear Sciences Limited, Lower Hutt, New Zealand and SeismicUnix from CWP (Center for the Wave phenomena, Colorado School of Mines). We thank all the people involved directly or indirectly in the elaboration of this work.