Characterizing a Naturally-Fractured Carbonate Formation for a CO 2 Storage Operation

Investigation into geological storage of CO2 is underway at the Technology Development Plant (TDP) at Hontomín (Burgos, Spain), the only current onshore injection site in the European Union. The storage reservoir is a deep saline aquifer located within Low Jurassic Formations (Lias and Dogger), formed by fractured carbonates with low matrix permeability. Understanding the processes involved in CO2 migration within this kind of low-primary permeability carbonates influenced by fractures and faults is key to ensure safe operation and reliable plume prediction. During the hydraulic characterization tests, 2300 tons of liquid CO2 and 14000 m of synthetic brine were co-injected on site in various sequences to characterize the pressure response of the seal-storage pair [de Dios et al, 2017] The injection tests were analyzed with a compositional dual media model which accounts for both temperature effects (as the CO2 is liquid at the bottom of the wellbore) and multiphase flow hysteresis (to effectively simulate the alternating brine and CO2 injection tests that were performed). The pressure and temperature responses of the storage formation were historymatched mainly through the petrophysical characteristics of the fracture network [Le Gallo et al, 2017]. The dynamic characterization of the fracture network dominates the CO2 migration while the matrix does not appear to significantly contribute to the storage capacity. This initial modeling approach was improved using the characterization workflow developed within the European FP7 CO2ReMove project for sandstone fractured reservoirs [Ringrose et al., 2011; Deflandre et al., 2011]. Fractured reservoirs are challenging to handle because of their high level of heterogeneity that conditions the reservoir behaviour during the injection. In particular, natural fractures have a significant impact on well performance [Ray et al, 2012]. Furthermore, the understanding of the processes involved in CO2 migration within relatively low-permeability storage influenced by fractures and faults is essential for enabling safe storage operation [Iding and Ringrose, 2010]. As part of the European H2020 ENOS project, the site geological model is updated by integration of the recently acquired data such as the image log interpretations from injection and observation wells. The geological model is generated through the analysis and integration of data including borehole images and well test data. Following a methodology developed for naturally fractured hydrocarbon reservoirs [Ray et al., 2012], the image log analysis identified two sets of diffuse fractures. A Discrete Fracture Network [Bourbiaux et al., 2005] was built around both wells which encompass the caprock, storage and underburden formations. The fracture characteristics of the two sets of diffuse fractures, such as orientations, densities and conductivities, are calibrated upon the interpretation of the injection tests [Le Gallo et al, 2017]. For each facies, the DFN characteristics were then upscaled and propagated to the full-field reservoir simulation model as 3D fracture properties (fracture porosity, fracture permeability and equivalent block size). Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 23 May 2018 doi:10.20944/preprints201805.0324.v1 © 2018 by the author(s). Distributed under a Creative Commons CC BY license. 2


Introduction
CO2 geological storage has reached industrial scale in sites such as Sleipner (Norway), In-Salah (Algeria), and Decatur (USA).These sites represent examples approaching the ideal conditions for establishing a commercial site according to criteria established in SACS project [Chadwick et al. 2007].Other CO2 injection pilots also achieved notable scientific results for example including Otway (Australia), Ketzin (Germany), Nagaoka (Japan), Lacq (France) [Kovács et al, 2015].The Hontomín pilot is the only current onshore injection site in Europe for CO2 geological storage, recognized by the European Parliament [EP resolution 2014] as key test facility for CCS technology development.It is located close to Burgos in the north of Spain, owned by Fundación Ciudad de la Energía (CIUDEN).The Hontomin storage reservoir comprises fractured limestones and dolomites .The injection of CO2 into fractured carbonate rocks for dedicated storage is unique in a European context, although considerable experience of injection into carbonates was gained in North America in association with enhanced oil recovery operations as for example in the Weyburn oil field [Whittaker et al, 2011].
Prior to dynamic modeling the CO2 migration, a detailed modeling of the characteristics of the storage complex is required as a key first step workflow for CO2 storage site characterization [Delprat-Jannaud et al., 2015].This paper presents the main features of the geological model with its stratigraphic and petrophysical properties.This work integrates the characterization (size, conductivities) of the natural fracture networks by its modeling in the vicinity of the wells [de Joussineau et al, 2016].

Geological context
The storage site represents a structural dome where the cap rock and reservoir belong to the Jurassic formations Marly Lias and Sopeña respectively.Keuper is the underlying seal and Dogger, Purbeck and Weald form the overburden [Rubio et al 2014].The representative lithological column of Hontomín site and geological cartography of the area are shown in the Figure 1.The Marly Lias and Pozazal formations are mainly comprised of highly carbonated marls (close to 50%) with uniaxial strength values equal or higher than 130 MPa and Young Modulus values in the range 15-30 GPa.The Sopeña formation is comprised of limestone in the upper part of reservoir and dolomite at the bottom.Their geomechanical properties correspond to rocks with high values of uniaxial strength, as the cap rock case, being the values for limestone and dolomite equal or higher than 180 MPa and 190 MPa respectively.As regards the Young modulus values, they are in the ranges of 30-60 GPa and 50-85 GPa [Kovács, 2014].These data suggest that there are different post fracture behaviors for the seal and reservoir, taking into account the tectonic effects induced in the rock layers during the formation of the dome.Uniaxial strength values are really high for a carbonated marl, but the Young modulus range reveals that despite being a strong rock its behavior should be of type "strain softening", with relevant deformations pre and post fracture.On other hand, carbonates (limestone and dolomite) show also high values of uniaxial strength, but Young modulus values support the behavior should be of type "elastic brittle" that means no significant strains take place after rock fracture.
According to mentioned assumptions, Hontomín geological model covers two different fault groups which affect the reservoir and overburden respectively, but there is no continuity between them through the cap rock what ensures its integrity, except two cases described below.Figure 2 shows the arrangement of faults and fractures in the study area.On the other hand, the values of some petrophysical properties determined by laboratory tests, as effective porosity, gas permeability and injections of CO2 and brine in reservoir conditions to assess the hydrodynamic and geochemical effects, revealed that fluid transmissivity is dominated by carbonate fractures, while the matrix does not appear to significantly contribute to the storage capacity.Results from subsequent tests conducted at field scale on site during the hydraulic characterization phase supported this assumption [de Dios et al, 2017].
Two main faults cross the storage complex from the reservoir to the overburden, which are considered limits in the geological model.Ubierna fault, located at the southern part of Hontomín area, and East fault are shown in Figure 3. Future works will be carried out to determine if both faults are sealing or transmissible to flow.

Geological model
The geological model covers the whole storage complex from the overburden (Dogger formation) down to the storage (Sopeña formation) and the underburden (Keuper formation) as shown in Figure 1 Figure 1.Vertical cross-section of the geological model of the Hontomìn storage complex; The HA and HI wells are the observation and CO 2 injection wells respectively while H2 to H4 wells are legacy wells

Structural model
A 3-D seismic reflection survey was acquired in 2010 which parameters included 22 source lines (crosslines), deployed E-W, perpendicular to 22 receiver lines (inlines) deployed N-S, with intervals of 25 m between sources and between receivers; the inline and crossline spacing was 250 m and 275 m, respectively, covering a total extent of 36 km 2 [Alcalde et al., 2014].Due to the complex geological setting of the Hontomìn site and the existence of an unexpected sharp velocity inversion near the surface, associated to the Upper-Lower Cretaceous contact [Alcalde et al., 2014], the 3-D seismic only identified the main horizons below the Dogger.These horizons were matched to the corresponding makers for legacy and newly-drilled wells.
The faults are only interpreted at the top of the storage formation from the 3-D seismic interpretations and are assumed to be vertical.
The grid of the geological model was designed to follow the facies vertical heterogeneities.
The lateral heterogeneities representation are coarser due to lack of well correlations.

Petrophysical model
The petrophysical model is established upon the facies and porosity log available from HA and HI wells and laboratory results.Since only the two newly-drilled wells have facies and porosity information, a simple modeling approach is selected for property modeling.

Facies
The vertical correlation length of the facies is assumed to be a function of the formation thicknesses.This is particularly important in the Pozazal, Marly Lias formations which show successions of shale, limestone and marls.The vertical correlation length may alter the vertical connectivity.Consequently, the grid thickness is quite small in the Pozazal formation which shows alternances of marl and shale (see Figure 4).

Porosity
The porosity is distributed within the grid based upon a moving average algorithm from the up-scaled porosity at the newly-drilled wells ensuring they are constrained to the facies distribution in the grid.

Permeability
The permeability is assumed to be constant per facies.The matrix permeability is mainly available from core measurements [Kovács, 2014;Dávila et al, 2016] which may be complemented for some facies by literature data [Bennion and Bachu, 2007] which is consistent with the core measurements.The previous dynamic modeling [Le Gallo et al, 2017] show a minor contribution of the matrix in the history matching of the pressure response.The fracture permeability is the main driver to the history matching of the pressure response [Le Gallo et al, 2017].

Fracture characterization and modeling
A study on fractures characteristics in both wells was conducted during the well logging campaigns, which particularly focused on the reservoir strata [de Dios et al, 2017].Data gained from acoustic televiewer during the performance of these works provided information about the type of fractures and their distribution [Collier and Ridder, 1992].

DFN modeling
As summarized by Jing [Jing, 2003], the DFN method is a special discrete model that considers fluid flow and transport processes in fractured rock masses through a system of connected fractures.DFN are particularly well suited the applications on characterization of the permeability of fractured rocks and generic studies of fracture influences, and the design of rock engineering works for near-field problems.
The stochastic simulation of fracture systems is the geometric basis of the DFN approach and plays a crucial role in the performance and reliability of the DFN model.The key process is to create PDFs of fracture parameters relating to the densities, orientations and sizes, based on field mapping results using borehole logging data and scanline or window mapping techniques, and generate the realizations of the fractures systems according to these PDFs and assumptions about fracture shape (circular discs, ellipses or polygons).The approach implemented in FracaFlow™ is based upon a pipe model which represents a fracture as a pipe of equivalent hydraulic conductivity starting at the disc center and ending at the intersections with other fractures, based on the fracture transmissivity, size and shape distributions [Cacas et al, 1990] To include the fractures identified in the televiewer log, a specialized software, FracaFlow™ [de Joussineau et al, 2016] is used to better integrate the different fractures into the flow model and calibrate their conductivities on the well test interpretations.Depending upon the geological context, different types of objects may be identified at different scales (Diffuse fractures, fault-associated fractures, faults) which may be intercepted by the wells [Ali et al, 2009].In general, horizontal wells provide the best source of information on the fractures.However, at Hontomìn, all wells including the newly-drilled one are sub-vertical wells.
The approach relies upon a DFN construction from the identified fracture and calibration of the various parameter upon the well tests [Ray et al, 2012].A fracture is characterized by its attributes such as orientation (dip and strike), length, aperture, and origin or morphology.A fracture set is a group of fractures having similar attributes.The fracture set can be described with the previous attributes and the density of fractures.The fracture characteristics may be associated to the facies generated on the grid.The grid and its properties (facies, porosity, and permeability) for the matrix, the main horizons (Top Dogger, Top Lias, Top Calcareous Sopeña, and Top Keuper), the fault maps at Top Calcareous Sopeña and the well trajectories and logs (facies, porosity) are imported into FracaFlow™.All the fractures characteristics for HI and HA from televiewer log [de Dios et al, 2017] were also imported as fracture logs (type, dip and azimuth as a function of depth).As shown in Figure 4, 5 types of fractures (0,10,14,16,18) were identified from the Dogger to the well shoes (Sopena/Carniolas formations).

Fracture analysis
The pole analysis is based the statistical analysis of the dip and strike of the fractures.The available fracture may be analyzed based upon the facies which does not allow clear discrimination between the different fracture families.The pole analysis (Figure 5) is based on the statistical analysis of the dip and strike of the mean pole of each fracture family.The mean poles are analyzed according to their main orientations (dip, strike) as shown in Table 1.Most of the fractures identified in the wells (Figure 5) represent the factures associated with type '0' which are sub-horizontal and reflects the formation bedding.Two main family of diffuse fractures are identified (Table 1): one with an approximate north/south orientation (strike ~176 N) as shown in Figure 6, one with an approximate East/West orientation (strike ~85N) as shown in Figure 7.They regroup the fractures identified on the televiewer analysis associated with 10, 14, 16, and 18 types from Figure 4.The diffuse fracture families are then associated with the facies in terms of fracture densities as shown in Table 2.There is no fracture identified in the anhydrite facies but all the other facies exhibit some level of facture.
The estimated thickness of the various beds may be obtained from the density of the bedding fracture of the reservoir which is about 1 m for the limestone and dolomite.The parameters for size distribution, aperture, and conductivity are fairly uncertain and poorly constraint at this stage.The conductivity will be matched on the available well test.The parameters considered in the automated KH calibration are the conductivity of E-W diffuse fracture, and conductivity of N-S diffuse fractures.In order for the calibration algorithm to converge a large enough number of initial solutions, i.e. greater than 500, must be computed.The results reported in Table 3 are computed for at least 1000 generations of the genetic algorithm which is significantly beyond the expected default convergence.The calibration is performed on Upper Calcareous Sopeña with the matrix properties described previously (see section Petrophysical model) As shown in Table 3, convergence of the algorithm was only obtained on the last of the well test for both diffuse fracture sets given the DFN generated, the random number used to generate the analytical calibration.This value will be used to validate the DFN properties based upon the Upper Calcareous Sopeña properties.The results obtained for the March 2015 test indicates that a horizontal anisotropy is generated by the diffuse fracture network: the horizontal anisotropy ratio is about 5 (Table 3) between North-South and East-West directions.
The automated KH calibration options of FracaFlow™ are designed for conventional fractured reservoirs.The computed conductivities of the fracture networks are very low and the analytical automated KH calibration does not converge.The computed conductivity of the DFN around both HI and HA wells are shown in Figure 9.
The final step is the upscaling of the DFN properties (fracture conductivities) over the whole model.Upscaling enables to calculate equivalent parameters linked to diffuse fractures (permeability, porosity, block dimensions) using an analytical upscaling method.Consequently, the fracture properties will be spread to all facies (except anhydrite) and all zones of the geological model.Equivalent properties are then be computed for the fracture media at the scale of the grid block while the matrix properties are not modified.The matrix properties are not altered.This model of the storage complex of the Hontomìn will later be used as a basis for the history matching of the dynamic model for the injections planned in ENOS project and in subsequent works to be conducted on site.

Figure 1 .
Figure 1.Lithological column and geological cartography of Hontomín area

Figure 2 .
Figure 2. Faults and fractures within reservoir and overburden at Hontomín site

Figure 3 .
Figure 3.South (Ubierna) and East faults are limits of Hontomín geological model

Figure 4
Figure 4 Facies, fracture orientation-televiewer log (Courtesy of ICTJA-CSIC)-and porosity logs for HI and HA

Figure 5
Figure 5 Poles stereo diagram for the bedding (blue), N-S (orange) and E-W (green) fracture families

Figure 6 N
Figure 6 N-S fracture characteristics: fracture density (left), strike (center) and dip (right) diagrams with confidence ellipses for the population mean

FaciesFigure 8
Figure 8 Facies, fracture dip (bedding in dark green), interpreted N-S (purple) and E-W (light green) fracture families and porosity logs for HA DFN modeling

Figure 9
Figure 9 Calibrated conductivity (mD.m) of the E-W (blue) and in N-S (red) diffuse fracture networks in DFN around HI well in yellow (top) and HA well in green (bottom) in the storage formation

Table 2
Mean fracture density (Number of fracture /meter) as a function of facies.

Table 3
Simulated conductivity in different fall-off for HI as reported by Le Gallo et al[2017].
(*) NC: no convergence for a minimum fracture conductivity of 10 -16 mD.m