Open Access
Numéro
Oil Gas Sci. Technol. – Rev. IFP Energies nouvelles
Volume 72, Numéro 4, July–August 2017
Numéro d'article 22
Nombre de pages 24
DOI https://doi.org/10.2516/ogst/2017021
Publié en ligne 6 septembre 2017

© A. Estublier et al., published by IFP Energies nouvelles, 2017

Licence Creative CommonsThis is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Introduction

Since 2001, French national authorities and public research institutions have been strongly committed in the field of CO2 geological storage, as a promising means in the struggle against climate change (e.g., Brosse et al., 2010a). In this context, the deep saline aquifers existing in the Paris Basin formed the main target of successive research studies, in particular the PICOREF and Géocarbone projects the results of which were published in two issues of Oil & Gas Science and Technology 2010 (No. 3 and No. 4).

At this time, a large area located in the southeast of Paris was the focus of attention, because it is the region of the basin where most of the geological information, acquired during exploration for petroleum or geothermal energy, was available. In this area, the Dogger aquifers, notably the Oolithe Blanche carbonate formation, were the main considered aquifers (Brosse et al., 2010b; Delmas et al., 2010).

Following these research studies, a 4-year demonstration project driven by industry (Bader et al., 2014) was initiated in 2008. It aimed to study the possibility to set up a pilot infrastructure of CO2 transport and storage in the Paris Basin. The first objective of the France Nord project was to identify a geological site in the deep saline aquifers of the Paris Basin providing a storage capacity of at least 200 Mt of CO2 over 50 years of injection. The project also aimed to improve and validate the methodology for selection, characterization and capacity estimation of storage sites for facilities of industrial size. The chosen site had to satisfy the following criteria:

  • the geographical and geological features of the site, and the petrophysical properties of the reservoir had to guarantee the safe storage of large volumes of CO2 over a very long time;

  • the reservoir had to exceed 1,000 m depth (below ground level) to be sure CO2 is injected as a supercritical fluid and storage volume is thus minimized;

  • the water of the aquifer had to be undrinkable and unsuitable for agricultural use due to its salinity (salinity > 10 g/L);

  • finally, any risk of interaction with potential surface or subsurface activity had to be avoided.

Only three areas of the Paris Basin located in the Triassic sequences (two areas, North and South, in the Keuper and one in the Buntsandstein) were finally selected for the characterization of the storage potential by optimizing the number of injection wells in order to inject 200 Mt of CO2 over 50 years.

Next to this industrial project, a R&D project funded by ADEME (French Environment and Energy Management Agency) and several industrial partners (TOTAL, ENGIE, EDF, Lafarge, Air Liquide, Vallourec) was initiated to analyze the risks associated to the selected pilot CO2 storage.

For that purpose, the R&D project was based on the European Directive (Directive 2009/31/EC on the Geological Storage of Carbon Dioxide) – hereafter the Directive – which defines the basic geological entities of a storage site, and stresses the importance of modeling for site characterization and risk assessment. According to Annex I of the Directive, a numerical simulation based on a “three-dimensional static geological earth model” is able to predict the “Storage dynamic behavior” (ibid.). The approach appears essential to establishing a monitoring plan and to evaluating risk scenarios, which are two mandatory aspects, but also to optimize injection operations and reservoir occupancy. Thus, for a selected area of the Paris Basin that could potentially hold storage site, a fluid-flow model associated to a geological model was built. Regarding the security issue, the injected CO2 in its integrity must remain confined in the geological portion of what is called the storage complex. This entity plays a central role in terms of security, as it represents the place from which CO2 should not escape, even dissolved in water. When hydrodynamic trapping drives the behavior of the CO2 plume, modeling is the only way to delimit the storage complex at an initial state of the project. Moreover, risk analysis study derives from a clear understanding of storage impact, which in turn is appraised through numerical modeling of multiple relevant phenomena. In the France Nord project, the risk assessment was focused on the impact of the pressure build-up at the basin scale and on the geochemical impacts in the storage complex.

This article presents the characterization of the selected area (South Keuper) studied in the industrial project and the associated risk assessment carried out in the R&D project. Thus, after presenting the geology of the Paris Basin and the potential of the CO2 storage site, the geological model and the fluid flow model are described. Then, an optimization of the number of wells and a sensitivity analysis to some uncertain parameters are presented and discussed. Then, the geochemical effects induced by CO2 injection in the Triassic formation are analyzed, both in the short (injection period) and long term (storage period of 1,000 years). This goal is achieved by simulating the fluid-rock interactions with a coupled geochemical and fluid flow model in the storage complex.

This work was performed using the research simulation software COORES coupled with ARXIM. COORES (Le Gallo et al., 2006) is a 3D multiphase code of IFP Energies nouvelles (IFPEN), which allows to simulate fluid flow and species transport. The model governing equations, solved simultaneously, are the mass conservation of each element, the Darcy’s law and the energy conservation. ARXIM (Moutte, 2009; Moutte et al., 2010) is a 0D geochemistry code that allow to model fluid-rock interactions. It allows to simulate over the time the evolution of mineral proportions and of water composition. It was developed by the Saint-Étienne School of Mines (http://www.emse.fr/~moutte/arxim/).

1 The Potential CO2 Storage Site in the West Paris Basin

The Paris Basin (Perrodon and Zabek, 1990; Guillocheau et al., 2000; Gonçalvès et al., 2004), the largest sedimentary basin in France, represents the French part of the intracratonic Anglo-Paris Basin. All its limits are erosional. A simple layer stacking model is not relevant in this kind of basin, where peripheral polygenic truncations due to frequent tectonics events. In such an epiorogenic basin the hydrogeological behavior should be determined by either cross-boundary effects or fault-leaking effects. The basin covers ca. 140,000 km2 onshore and ca. 10,000 km2 offshore in the French sector of the Channel.

About one third of the 6,000 exploration wells drilled in the Paris Basin reached either Paleozoic series or the base of Jurassic ones. A reliable post-Permian geological sketch has been drawn from this database (Perrodon and Zabek, 1991; Guillocheau et al., 2000). The post-Paleozoic sedimentary pile reaches its maximum thickness, ca. 3,500 meters, in the Brie east area of Paris.

The main tectonic and depositional events of the basin occurred in the selected area. A comprehensive stratigraphic column is provided by the lithological log of the Tousson 101 (TSN 101) well, located just in the middle of the studied sector (Fig. 1). The well drilled through the Mesozoic and the thin outcropping Tertiary sequences. Jurassic and Cretaceous deposits have the highest thicknesses and consist mainly of carbonate sediments; two siliciclastic intervals are representative of the Lower Cretaceous and the Triassic sediments (the target stratigraphic interval). No long lasting hiatus interrupts the sedimentation in this rather preserved area.

thumbnail Figure 1

Stratigraphic column of the South Keuper sector area.

Several main tectonic faults are recognizable on the global structural map of the top of the Triassic (Fig. 2, Delmas et al., 2002). They could represent rejuvenated tardi-hercynian faults, as Pays de Bray or Étampes ones, or even younger faults, as the N40 ones related to the Liassic distension. The map of Figure 2 corresponds to the top of the Triassic geological model described further. An orange meridian line marks its eastern limit.

thumbnail Figure 2

Structural map at the top of the Triassic formations. The potential CO2 storage sites should be located in the areas represented by either the black rectangle (South Keuper) or the yellow one (North Keuper). The orange meridian line is the eastern limit of the geological model.

The best choice for a massive CO2 geological storage consists of the Upper Triassic deposits situated 70 km southwards of Paris, in a perimeter called “South Keuper sector” (Fig. 2 and 3). In this area, the Upper Triassic or Middle Keuper series are characterized by continental redbeds deposited in subtropical to subarid climates. It belongs to the western part of the Germanic facies realm (Bourquin et al., 2002).

thumbnail Figure 3

Hydraulic layers. The South Keuper sector area is located in the yellow rectangle.

The selected sector (framed in black in Fig. 2), located near a depocenter of secondary importance, called the Pithiviers graben, contains two types of faults. The first one, illustrated by the Étampes fault, affects all the Mesozoic and Eocene series, while the other type of faults, oriented N40, is sealed at the top of the Middle Jurassic. In the sector, the tectonic style is controlled by a N40 parallel stripe fault system, which interferes westward with a deep sub meridian graben along the Étampes fault. Down to the Triassic levels, the structural traps are typical of a distensive context, expressed as tilted fault-bounded closures.

1.1 Characteristics of the Storage Complex

Based on their hydraulic behavior at the regional scale (aquifer vs. aquitard vs. aquiclud), the Triassic beds were subdivided in eight “hydraulic layers” according to a previous modeling study of the regional aquifer (Houel and Euzen, 1999) (Fig. 3).

Beds list and names from the bottom to the top of the Triassic serie:

  • cap rock: the Marine Rhaetian,

  • the Upper Boissy sandstones/Upper Châlain marls,

  • the Lower Boissy sandstones/Middle Châlain marls,

  • the Lower Châlain marls,

  • reservoir (CO2 storage): the Chaunoy sandstones/Upper variegated marls,

  • the Intermediate Sandstones and Clays,

  • reservoir (CO2 storage): the Donnemarie Sandstones,

  • the Bundsandstein Sandstones.

The immediate availability of a complete data set directly adapted to a first assessment of the CO2 migration explains the choice of this relatively low-resolution layering.

Among these eight layers, the fifth one called “Chaunoy sandstones-Upper variegated marls” was selected as the best candidate to massively store CO2 in subsurface. In the South Keuper sector, it is a 60 m-thick, multistorey, siliciclastic reservoir deposited along a piedmont alluvial apron passing eastward to an evaporitic continental mud-flat or bajada (Spötl and Wright, 1992). As shown in a detailed section of the Middle Keuper of the TSN 101 well (Fig. 4), the reservoir levels are composed of dolomite-cemented subarkosic sandstones. They are interbedded with anhydritic red marls. Thick massive to vuggy dolomite, related to groundwater dolocretes, often caps this large sand sheet. The alluvial apron is 50 km wide. The sandstone net-to-gross ratio decreases eastwards. It is pretty high in the sector center (0.88), even with a porosity cut-off of 5% (0.63), but the mean porosity is quite low with values around 10%, probably explained by deep burial and the diagenetic history linked to highly saline brines (80 g/L NaCl). In this area less buried reservoirs of large extension contain potential water resources.

thumbnail Figure 4

“Quick-Look” layering of the Norian Chaunoy sandstones formation.

The subarkosic nature of the detrital content, associated with a complex diagenesis, explains the great diversity of the mineralogy of the Chaunoy sandstones. This variability appears clearly with the Quick Look method (Serra, 1985) applied to the TSN 101 well (Fig. 4), where the repetitive interlayering of sandstone, dolomite and marl beds is observed. The lateral extension of the drains is probably hectometric to kilometric. The dominance of sheet-like bodies, typical of the alluvial apron architecture, can significantly reduce the vertical connectivity. Only moderate faulting and local channeling can interrupt the continuity of the extensive clays seals or the widespread cemented paleosoils. The global architecture of the deposits is very similar to the one explored in the Chaunoy oilfield, located 30 km north-westward in the same apron (Comité des techniciens ELF, 1991; Eschard et al., 1998). In this large oilfield, the Chaunoy sandstones constitute the oil reservoir, trapped in a horst structure where production data revealed an unique Oil Water Contact (OWC). The structural closure, ca. 60 m, is equivalent to the formation thickness. We assumed here that the Chaunoy sandstone formation behaves as a single hydraulic unit in the studied area.

The Rhaetian dolomitic Châlain red marls form the immediate cover of the Chaunoy sandstones. They are cut by the fluviatile meandering Boissy sandstones, occurring as sporadic good reservoir levels. In the Lower Lias cover only the Hettangian series were integrated in the geological model as a hydraulic level. The ultimate cover is represented by the 200 m-thick section of the Lower Liassic marine marls. The integrity of this cover is proved by the trapping of oil in the Triassic reservoirs.

2 Description of the Models

Three models were used in this study. A large part of the Paris Basin was included in both the geological and the fluid-flow models because: (1) the hydraulic circulation of the Paris Basin is complex due to multiple faults, outcrops and an hydraulic connection between the Triassic and the Channel, (2) one of the study objectives is to simulate the overpressure impact at the basin scale since it may extend several tens or hundreds kilometers. In contrast, the geochemical model is limited to the storage complex. A one-way coupling method, explained in Section 4.1, was used to transfer at this scale the results obtained from the basin-scale models.

2.1 Geological Model

The Triassic geological model (Lecomte et al., 2010) used in this study covers a large part of the Paris Basin – approximately 400 km-500 km. It includes the eight geological Triassic units described in the previous section; among them, the Upper Boissy, the Lower Boissy, the Chaunoy and the Donnemarie formations are the four main reservoir units. The mesh contains 176,352 gridblocks (including 15,621 active cells) of 3 km-side length (132 in X-direction, 167 in Y-direction and 8 in Z-direction) (Fig. 5).

thumbnail Figure 5

Mesh of the geological model and schematic view of the boundary conditions, for the selected area: imposed pressure at the outcrops (0.1 MPa at 200 m subsea, in green, and 0.1 MPa at 0 m, in purple). The mesh represents the top of the Triassic formations. It is limited at the eastern boundary by the meridian line defined in Figure 2.

The Paris Basin is structured by multiple faults (Fig. 6a). Only faults (Fig. 6b) whose throw exceeds 50 m in the Triassic were considered in the geological model. However, the hydrodynamic impact of faults with a small throw could be crucial in a small-scale approximately modeling, particularly around an injection site. In a basin-scale hydrogeological model it seemed a reasonable approximation to represent a fault by its throw value and to neglect the fault heave.

thumbnail Figure 6

Structural map showing faults and the Triassic top layer: a) faults (black and pink lines) of the Paris Basin; b) faults of the geological model illustrated in Figure 5.

The 3D faulted model of the Triassic of the western Paris Basin has been digitized from the Regional Atlas of Western Paris Basin (Houel and Euzen, 1999). Approximately forty maps, which include data attached to each of the eight hydraulic units (total thickness, effective thickness, porosity, permeability, and salinity), were used in addition to the previous faults pattern.

Porosity maps were obtained from the interpretation of well logs and the regional interpolation was based on geological interpretations. Permeability maps were computed using permeability/porosity correlation laws coming from petrophysical analysis of all cored wells of the western Paris Basin. The temperature values at the top of the Triassic sediments were provided by a BRGM database (BRGM, 2008) and a temperature gradient of 0.0345 °C/m was applied to compute the temperature in each layer.

Pressure maps from the Regional Atlas of Western Paris Basin (Houel and Euzen, 1999) are iso-pressures. Values of piezometric maps (relative to the sea level) were deduced from iso-pressure maps taking into account the variable salinity of Triassic waters.

Triassic reservoirs in Western Paris Basin do not represent a closed system. In the south part of the basin (between the towns Saint-Amand-Montrond and Montluçon), outcrops of the Triassic at about 200 m above the sea level over an east-west oriented area (about 130 km long and 6 km wide, Fig. 7) are considered as an area of hydraulic recharge at an average altitude of 200 m (Boissy and Chaunoy formations). Moreover, a hydraulic connection between the Triassic and the Channel was assumed in this work (from Boissy to Donnemarie formations). Indeed, a SE-NW gradient (towards Dieppe) of salinity in the Chaunoy layer leads to salinity values close to that of seawater.

thumbnail Figure 7

Triassic outcrops (purple area) at the west of the Paris Basin. Overlay of cultural data, surface topography (radar image) and isobaths (m) at the top of the Triassic series.

Table 1 summarizes the mean value of the properties of the geological model (Fig. 8) for each formation modeled in this study. As the Marine Rhaetian, the Lower Châlain marls and the Buntsandstein sandstones have a small thickness, they were defined as dead cells (porosity equal to zero and horizontal permeability equal to zero) in the geological model. For the vertical permeability, the 0.1 ratio, which is commonly used for sandstones, was assumed.

thumbnail Figure 8

Structural grid with porosity values and location of the injection well defined in Section 3.1.

Table 1

Mean value of the properties for each geological formation modeled for this study.

2.2 Fluid Flow Model

2.2.1 Petrophysical Properties

In addition to porosity, horizontal permeability, net-to-gross, salinity and temperature values provided by the geological modeling, the fluid flow model required rock compressibility, vertical permeability, relative capillary and capillary pressure properties.

Average rock compressibility was estimated using existing values for sandstones (Newman, 1973). It is 4 × 10−10 Pa−1. The fluid relative permeability (K rw and K rG) and capillary pressure P c (Fig. 9) are computed with the Brooks-Corey local model (Brooks and Corey, 1966). Correlations were used to determine the two parameters of the model (Fig. 9) using an average permeability value:

  • the local displacement pressure, noted $ {P\mathrm{c}}_{\mathrm{d}}$ (threshold pressure of the capillary pressure in drainage mode) (Thomas et al., 1968);

  • the parameter λ related to the index of pore size distribution.

thumbnail Figure 9

Water and gas a) relative permeability curves and b) capillary pressure (0.1 MPa) curve.

  ( P C d P c ) λ = S W - S Wi 1 - S Wi = S W * ;   K r W = ( S W * ) 2 + 3 λ λ ; K r G = ( 1 - S W * ) 2 ( 1 - ( S W * ) 2 + λ λ   ) $$ {\enspace \left(\frac{{P}_{{\mathrm{C}}_{\mathrm{d}}}}{{P}_{\mathrm{c}}}\right)}^{\lambda }=\frac{{S}_{\mathrm{W}}-{S}_{\mathrm{Wi}}}{1-{S}_{\mathrm{Wi}}}={S}_{\mathrm{W}}^{*};\enspace {K}_{r\mathrm{W}}={\left({S}_{\mathrm{W}}^{*}\right)}^{\frac{2+3\lambda }{\lambda }};\hspace{1em}{K}_{r\mathrm{G}}={\left({1-S}_{\mathrm{W}}^{*}\right)}^2\left(1-{\left({S}_{\mathrm{W}}^{*}\right)}^{\frac{2+\lambda }{\lambda }}\enspace \right) $$with the water saturation S w and the irreducible water saturation S wi (= 30%).

2.2.2 Thermodynamic Model

Pure CO2 is assumed to be injected in a saline aquifer with a variable salinity. Thermodynamic properties of the injected fluid are listed in Table 2.

Table 2

CO2 thermodynamic properties (Reid et al., 1987).

Water density and viscosity are computed from correlation models (Schmidt, 1969) based on pressure, temperature, salinity and the molar fraction of dissolved CO2. The Peng-Robinson equation of state (Peng-Robinson, 1976) calculates the gas density and the Lorenz-Clark equation (Lohrenz et al., 1964) estimates the gas viscosity. BOth uses the CO2 thermodynamic properties defined in Table 2

CO2 is injected into the Chaunoy formation in the South Keuper. At the injection depth, i.e., 2,300 m depth, temperature is close to 93 °C (366.15 °K) and pressure about 25.5 MPa. At these temperature and pressure conditions, CO2 is a supercritical fluid. Its viscosity (4.7 × 10−15 Pa s−1) is close to the viscosity of the gas phase and its density (608 kg/m3) is close to the density of the liquid phase.

CO2 dissolution is modeled only when fluid-rock interactions are simulated. The CO2 dissolution in brine is a relevant effect for pressure decrease during CO2 injection and thereafter; not considering this effect introduces an error in the simulations leading to an overestimation of the pore pressure. However, the discretization in space of the fluid flow equations leads to numerical diffusion in case of dissolution modeling; this causes an overestimation of the amount of dissolved CO2 in the first years of injection and therefore, an underestimation of the pore pressure. The mesh refinement leads to the decrease of the numerical diffusion. However, the fluid flow model was primarily built for a first assessment of the number of wells needed. CO2 dissolution has been considered in the subsequent analysis applied on a reduced volume on the geological model with a refined mesh. For this purpose, CO2 dissolution has been considered in the subsequent analysis applied on a reduced area of the geological model and with a refined mesh. For this purpose a sensitivity analysis to the CO2 dissolution was carried out and reported in Section 3 to qualitatively analyze the impact of its uncertainty on the injectable CO2 mass. In this case and for fluid-rock interaction simulations, CO2 dissolution is modeled by a table of equilibrium constants at different temperature and pressure. These constants are computed by the model of Soreide and Whitson (1992).

2.2.3 Boundary Conditions

In this study, two types of boundary conditions were defined:

  • imposed pressure (open boundary condition, Fig. 5) at the outcrops and at the potential connection with the Channel;

  • no-flow conditions on all other borders of the model.

Imposed boundary conditions at the outcrops (0.1 MPa at 200 m of altitude, in green in Fig. 5) imply an input of fresh water in the border cells. At the other boundary of the model, the existence of a connection with the Channel (condition: 0.1 MPa at 0 m, in purple in Fig. 5) is only supposed but it appears relevant to create a permanent flow from the South to the North. However it is likely that this flow of water does not only come from the Triassic formations but also from the upper aquifers that are not modeled in our case.

An overpressure value equal to half the initial pressure value was considered, according to a commonly-used rule from the gas-storage industry. The fluid rate corresponding to this maximum pressure at the uppermost perforation of the injection well(s) leads to determine the number of wells which are needed to inject 200 Mt of CO2 in 50 years.

2.2.4 Initial State

All the geological formations of the model are initially saturated with brine: oil and gas zones are not modeled in the present study. In most fluid flow simulators, the calculation of the initial state of pressure consists in a single hydrostatic gradient. This is done by ranking all the gridblocks of the model on a 1D column, by averaging temperature and salinity for a given depth and then by extending the obtained information on the 3D model. Such an approach is correct as long as there is no lateral variation of the brine density. However, in our Triassic model, there are lateral variations of the temperature and of the salinity that intrinsically implies lateral variations of the brine density. Therefore, it appears necessary, before any injection simulation, to first perform balancing of pressures. To consider both boundary conditions and lateral variations of brine density, a prior fluid flow simulation is carried out over a period of 10 million years without any injection until obtaining a permanent flow regime.

To validate the steady state conditions obtained at the end of the simulation, a comparison with the pressure field of the Paris Basin was carried out. Figure 10 shows both the simulated pressure values with those measured (189 values) on a large number (132) of monitoring wells. The good correlation between the simulated and measured data allows one to conclude that the simulated pressure field is consistent with that of the Paris Basin.

thumbnail Figure 10

a) Location of the control wells chosen in the Paris Basin in order to b) compare measured pressure values, P gauge, and simulated ones, P sim (right).

2.3 Geochemical Model

Based on a selection of available data on mineral paragenesis and water composition, the first step was to define a geochemical model, i.e., a representation of the mineral composition either present in the virgin reservoir sediment or resulting from interaction with dissolved CO2.

The geochemical modeling was limited to the Chaunoy sandstone because it is the only reservoir layer for which sufficiently chemical and mineralogical data are available. The chemical composition of pore water from Triassic sandstones is provided by the Final Reports publicly available for several wells in the site area (Beaumont N.101, Chaunoy N.101, Chaunoy N.102, Étampes N.1, Les Roches Moreau N.1, Tousson N.101, Videlles N.1, and Villemer N.102). Mineral composition of the Chaunoy Sandstone, based on petrographic observations, is available from Bennehard’s (1970) PhD, in particular from the well Villemer N.122 (vol.%):

  • 80% quartz, 15% alkali feldspar, 5% plagioclase in the coarse fraction,

  • 90% illite or muscovite, 10% illite-smectite mixed layers in the fine fraction,

  • carbonate fraction exclusively composed of dolomite,

  • presence of kaolinite-rich interlayers.

This mineral description was used to infer a most likely composition of natural water at reservoir conditions. Geochemical modeling (speciation, mineral-water interaction) was achieved with Arxim using the extended Debye-Hückel formalism for the activity coefficients. The reconstruction of in situ water composition was mainly based on the chemical data available from the well Beaumont N.101 (test N.2, 14th Feb. 1967, Triassic layers ca. 2,240 m logging depth – close to true vertical depth). The analysis was limited to the major elements, as it is common for this type of data (Tab. 3), thus it was under-constrained, particularly regarding the carbonate system.

Table 3

Water analysis provided for test N.2 at well Beaumont N.101 (Triassic layers ca. 2,240 m logging depth – close to true vertical depth).

Three steps were followed:

  • first step, at surface conditions (20 °C, 1 bar): a pH value of 6.6 is consistent with the measured alkalinity, and with a Total Dissolved Carbon (TDC) value of 0.017 mol/kg H2O;

  • second one, at reservoir conditions (82.4 °C, 250 bar): screening pCO2 values expected at this temperature (Coudrain-Ribstein et al., 1998) leads to a pH value of 5.47 and a TDC value of 0.0603 mol/kg H2O at equilibrium with disordered dolomite;

  • third one, discussing saturation with respect to alumino-silicates, at reservoir conditions: constraining dissolved silica by equilibrium with K-feldspar, and dissolved aluminium by equilibrium with kaolinite, the measured potassium value leads to an acceptable saturation state with respect to the minerals observed by Bennehard (1970). A TDC value of 0.0573 mol/kg H2O gives a pH value of 5.43 and respects equilibrium with disordered dolomite and calcite. Speciation also shows that this water is equilibrated with anhydrite and chalcedony, oversaturated (0.60 in log units) with respect to muscovite and undersaturated (−0.73) with respect to albite (“low-albite”). A moderate oversaturation with respect to quartz (here 0.23) is expected in formation waters at comparable temperature (e.g., Bazin et al., 1997; Kharaka and Hanor, 2003).

At this step, a dynamic calculation towards equilibrium with albite was tested, using kinetic rates provided by Palandri and Kharaka (2004). It leads to a heavily different water composition, which must be considered unlikely in the natural conditions of the reservoir.

The water composition obtained at third step was then charged with CO2, in order to mimic a natural reservoir water equilibrated with a CO2 fluid phase at 250 bar. The resulting TDC value is close to 2.0 mol/kg H2O, at pH 3.79. This water remains almost equilibrated with chalcedony, K-feldspar, kaolinite, disordered dolomite, calcite and anhydrite. At equilibrium with these minerals, the calculated values are the following: pH 4.62, TDC 2.11 mol/kg H2O, ionic strength 1.08, significant oversaturation with respect to muscovite (0.62) and undersaturation with respect to low-albite (−1.50).

Dynamic calculations involving either muscovite or albite lead to a profound alteration of the water-mineral system. Moreover, they reveal that the adopted mineral representation is inadequate to simulate a realistic evolution of the system. This difficulty was interpreted as a consequence of the too poor representation either of the mica/illite and plagioclase phases, likely to dissolve, or of the Na-rich clay minerals, likely to precipitate. In particular, the model does not consider solid solutions – unfortunately not available in the version of Arxim.

A remedy to this situation could be found using dawsonite as a sink for sodium. This hydrated carbonate mineral of sodium and aluminum (NaAlCO3(OH)2) has become popular in the context of geochemical studies dedicated to CO2 geological storage (Johnson et al., 2004). It has been observed in the sedimentary context (reviews in Baker et al., 1995; Hellevang et al., 2011; see also Gao et al., 2009), nonetheless quite rarely. Hellevang et al. (2011) and Kaszuba et al. (2011) investigated why geochemical modeling has a propensity to predict dawsonite formation, at least as a transient phase (e.g., Okuyama et al., 2013), whereas the mineral is rather scarce in nature or in experimental contexts. Activity diagrams show that an increase of pCO2 favors dawsonite stability, but can be counterbalanced by the decrease of pH value, which has the opposite effect. Dawsonite presents also a decreasing stability at increasing conditions of burial depth (increasing temperature and pressure). Hellevang et al. (2011) inferred from such considerations that in the temperature window which is the more favorable for dawsonite formation, kinetics has a large energy barrier to overcome, in other words precipitation requires a very high oversaturation level. On their side, Kaszuba et al. (2011) put forward that thermodynamic databases lack internal consistency between the mineral phases that potentially play a role in the relevant system (SiO2-Al2O3-CaO-Na2O-CO2), among which solid solutions should be added, in particular for smectites. We are perfectly aware of these limitations, but it was beyond the scope of the modeling work presented here to develop such a consistent and more comprehensive database. Thus, in addition to calcite, the geochemical model considers dawsonite as secondary minerals. The thermodynamic data taken into account for dawsonite were chosen consistently with the results published by Bénézeth et al. (2007). The precipitation rate of dawsonite is spoilt by a complete uncertainty. The same as for kaolinite is chosen, that can be considered as median between the high rates associated to carbonates, and very sluggish ones that characterize micas or quartz.

The kinetic rate law for dissolution and precipitation of the geochemical simulator Arxim is simplified from Lasaga (1984, 1995, 1998). It depends on the effective reactive surface area, the activation energy (E a) and the equilibrium state (p k ). Mineral grain is modeled by a floating sphere and the dissolution and precipitation phenomena are modeled by the sphere radius evolution. For the minerals which are initially missing, the initial radius is assumed to be 10−7 m. The composition and kinetic parameters considered for representing the reservoir in the geochemical simulations are given in Table 4.

Table 4

Composition and kinetic parameters considered for representing the reservoir in the geochemical simulations.

3 Storage Complex

In parallel with the technical projects, French public authorities carried out the necessary efforts to transpose into national law the European Directive on CO2 geological storage (EC Directive 2009/31, 23 April 2009). The transposition took place through a series of milestones, detailed in Perrette (2011) and was effective on October 31st 2011. In conjunction with the transposition, a working group was financed by the General Directorate of Risk Prevention, at the Ecology, Energy, Sustainable Development, Transport & Housing Ministry (MEEDDTL), with the objective of producing a technical guidance document framing the security of storage sites. This document (Bouc et al., 2011) specifies the vocabulary introduced by the Directive for characterizing a storage site, particularly in the case of the so-called hydrodynamic trapping – or Migration Assisted Storage (MAS) – which occurs when CO2 is injected in an aquifer which does not offer large volumes of structural trapping.

According to this document, the injected CO2 in its integrity must remain confined in the geological portion of what is called the storage complex. This entity plays a central role in terms of security, as it represents the place from which CO2 should not escape, even dissolved in water. According to the definition of the storage complex, it contains the reservoir, all the intermediary formations and the cap rock.

Here the storage complex was defined from a fluid-flow simulation assuming the best scenario in terms of injection strategy and injectivity of the wells. This simulation describes the extension of the reservoir where migration of the CO2 plume should occur, and the resulting overpressure field.

3.1 Fluid-Flow Simulation

The location of the injection wells has to satisfy two main requirements: (1) 200 Mt of CO2 must be injected in 50 years; (2) the injection rate is determined by a maximum pressure criterion detailed in Section 2.2. In addition, several other choices were formulated: the injection wells are vertical, have a 18 cm diameter, and drill the whole Triassic series. In the case where multiple wells are required, the maximum distance between wells is 30 km.

Two areas were initially selected for a potential injection into the selected Triassic formations, named South Keuper and North Keuper (Fig. 2) based on the following selection criteria:

  • the thickness of the sandstones,

  • no major fault affect the selected areas,

  • rather high value of water salinity.

For each area, a study domain was defined: a square of 21 cells of side (about 3,969 km2). In order to overcome any problem of direct interaction with faults, the cells that were too close to (less than 6 km) faults were excluded. A well was defined in one cell out of two. For each well (one by one), the storage capacity is tested according to the overpressure criterion.

The simulation work showed that the southern area is much more favorable than the northern zone to inject 200 Mt of CO2 over 50 years. Four wells would be sufficient in the south against twenty in the north because in the southern area the available reservoirs are the Chaunoy and Donnemarie formations whereas in the northern sector only the Chaunoy formation is present. These results were obtained using a base scenario, which is presented below.

The characteristics of the base scenario are listed below:

  • boundary conditions consist of constant pressure at the outcrops and at the Channel and no flux elsewhere;

  • vertical permeability anisotropy ratio equal to 0.1;

  • faults not sealed (multiplying transversal transmissivity factor equal to 1, longitudinal flow is not modeled);

  • no CO2 dissolution;

  • injection condition: imposed pressure equal to 150% of the initial pressure at the upper most sector of the four injection wells.

With these injection conditions, 264 Mt of CO2 is injected over 50 years that is beyond the scope (Fig. 11). The required amount of 200 Mt is injected in only 29 years. Among the four wells, Well_4 has the best injectivity: more than 100 Mt is injected over 50 years.

thumbnail Figure 11

Mass of injected CO2.

The CO2 injection in the Chaunoy formation causes an overpressure which extends over time (Fig. 12). After 50 years of injection, the observed 11.4 MPa corresponds to approximately 46% of the initial pressure value of Well_1. At the end of the injection, the overpressure above 0.1 MPa covers an area of 16,415 km2 (Fig. 13a) and the one above 1 MPa an area of 8,690 km2; 1,000 years after the end of the injection, the overpressure above 0.1 MPa covers an area of 51,869 km2 (Fig. 13b) and the one above 1 MPa an area of 28,810 km2. Therefore on the long term, the CO2 injection impacts the whole basin and it produces a flow at the outcrops and in the Channel due to hydrodynamism of the basin. The maximal overpressure (Fig. 13) decreases from 11.4 MPa at the end of the injection to 2.5 MPa at 1,050 years while the location of the maximal value moves towards areas where mobility is reduced, i.e., in the northern area of Chaunoy. The compartmentalization of the overpressure field shows that this one is strongly dependent on the network composed of faults with strong throw. However, to have a better description of the impact of faults on fluid flow, a more detailed model would be needed integrating all local (and not only regional) faults and defining faults as volume objects to simulate flow both across and along the faults. The CO2 plume extension (Fig. 14) remains low compared to the overpressured area: an area of 316 km2, at the end of the injection and of 524 km2, 1,000 years after.

thumbnail Figure 12

Maximal overpressure (MPa) in the Chaunoy formation over time (year).

thumbnail Figure 13

Overpressure (MPa) above 0.1 MPa at the a) end of the injection and b) 1,000 years after – Chaunoy layer.

thumbnail Figure 14

CO2 plume (gas saturation) at the a) end of the injection and b) 1,000 years after – Chaunoy layer. The limits of the storage complex are represented by the red rectangle.

The impact of the CO2 injection in the Keuper area on the exploited oil and gas fields was analyzed in terms of local overpressure. The simulated results are summarized in Table 5.

Table 5

Impact of the injection on exploited oil and gas fields.

Regarding the gas storages, no impact is visible because they are often outside the area impacted by the injection. Increase of the pressure of the oil fields would improve oil production conditions. The induced overpressure ranges from 0.9 MPa to 7.5 MPa.

In parallel of the risk analysis presented in this work, a second characterization study (Bader et al. 2014) was carried out after the refinement of the petrophysical properties maps and the integration of the secondary fault network. Finally, the grid definition was 1,000 m × 1,000 m with 28 stratigraphic layers (to be compared with the initial model 3,000 m × 3,000 m × 8 layers). This additional study revealed that the South Keuper area could hold the largest mass of CO2, nevertheless, at a too high cost in terms of injection-well number. In this article, only the simulation results before the model refinement are presented.

3.2 Sensitivity Analysis

A sensitivity analysis was carried out concerning the injection period in order to qualitatively analyze the impact of the uncertainty of the following input parameters on the injectable CO2 mass (Fornel et al., 2011, 2012):

  • maximum overpressure at injection point (uncertainty of 1% of the initial pressure);

  • porosity (uncertainty of 50% of the initial values);

  • horizontal permeability (uncertainty of 80% of the initial value);

  • rock compressibility (uncertainty of 50% of the initial value);

  • faults sealing (multiplying transversal transmissibility factor equal to zero or 1);

  • CO2 dissolution (modeled or not);

  • irreducible water saturation (uncertainty of 50% of the initial value) with dissolution modeling;

  • irreducible water saturation (uncertainty of 50% of the initial value) with no dissolution modeling;

  • vertical permeability anisotropy ratio (0.1 or 0.0001).

Each uncertainty range was arbitrarily defined based on expertise in the field of reservoir simulation. The initial value corresponds to the base scenario one. Each parameter has been characterized independently, without considering any interaction with the other parameters except for the irreducible water saturation and CO2 dissolution.

The results are presented in Figure 15 in a synthetic way in the form of two graphs: the relative sensitivity of a parameter and the relative error induced by a parameter on the mass of injectable CO2. The relative sensitivity of a parameter corresponds to the variation of the injectable CO2 mass relative to the base case, caused by a 10% variation of this parameter. It can be compared with the relative sensitivity of other parameters because it does not depend on the parameter uncertainty range. On the contrary, the relative error induced by a parameter on the mass of injectable CO2 corresponds to the variation of the injectable CO2 mass relative to the base case, caused by a variation of this parameter characterized by its uncertainty range.

thumbnail Figure 15

Sensitivities (red) and induced errors (green) of uncertain parameters on the mass of injectable CO2. (1) Maximum overpressure at injection point (uncertainty of 1% of the initial pressure), (2) porosity (uncertainty of 50% of the initial values, (3) horizontal permeability (uncertainty of 80% of the initial values), (4) rock compressibility (uncertainty of 50% of the initial value), (5) faults sealing (multiplying transversal transmissibility factor equal to zero or 1), (6) CO2 dissolution (modeled or not), (7) irreducible water saturation (uncertainty of 50% of the initial value) with dissolution modeling, (8) vertical permeability anisotropy ratio (0.1 or 0.0001), (9) irreducible water saturation (uncertainty of 50% of the initial value) with no dissolution modeling.

The comparison between the sensitivity graph and the error graph shows that the most sensitive parameter is not necessarily the parameter which induces the maximum of error and conversely for the less sensitive parameter because this depends on its uncertainty range. The most influential parameter is the maximum overpressure at injection point but it induces a small error because its uncertainty caused by the bottom pressure measurement is assumed to be only 1%. Similarly, the three parameters inducing the greatest errors are: horizontal permeability, porosity and fault sealing (in a decreasing order); whereas, the three most influential parameters are maximum overpressure at injection point, porosity, horizontal permeability (also in a decreasing order). The vertical permeability has a low impact and is not much sensitive. This is not a result that can be generalized because it is due here to a small mobility of the CO2 plume caused by the coarse mesh. In the case of an injection at the bottom of an aquifer with a dip or at the bottom of a system of communicating aquifers, results would be very different. Similarly, the irreducible water saturation and the CO2 dissolution are not much influential and cause a small error.

Whatever the uncertainty range value of all the studied parameters (except for two), the simulation succeeds in injecting 200 Mt over 50 years. Only the uncertainty of porosity and horizontal permeability prevents to satisfy the injection criterion in all cases. If the porosity or the horizontal permeability is reduced of more than 40% of the initial value, the injectivity is not sufficient.

For this preliminary site characterization, all uncertainties have not been studied such as the ones related to the geological model (geometry, gridding, secondary faults) and to the fluid flow model (relative permeability, hysteresis phenomenon, thermodynamics, fault modeling). Their no-negligible impact should be studied in a more complete site characterization.

3.3 Case of the France Nord Project

At this step of the study, the horizontal CO2 storage complex limits can be defined using the extensions of the CO2 plume (Fig. 15). The storage complex limits are assessed assuming that the stored CO2 must never cross them. The chosen complex storage area is four times bigger (two in each direction) than the simulated CO2 plume extension, assuming that the plume would be more extended at the top of the Chaunoy formation in case of a mesh vertically refined. A domain of 57 km × 69 km is defined. Vertically, the storage complex contains (1) the reservoir composed of all the sandstone formations, (2) the intermediary formations composed of the sandstone/marl formations and of the Lower Châlain marls and (3) the cap rock composed of the Marine Rhaetian. However, for a primary assessment of the fluid-rock interactions induced by the CO2 injection, only the geochemical effects in the sandstone formations will be simulated in the next section.

4 Geochemical Effects of CO2 Injection

Numerous studies simulating the induced fluid-rock interactions with a reactive transport simulator have already been carried out for other CO2 storage sites, for example: Dogger and Callovo-Oxfordian formations of the Paris Basin (Wertz et al., 2012), Utsira aquifer at Sleipner (Audigane et al., 2007), Upper Triassic Stuttgart Formation at Ketzin (De Lucia et al., 2014). However, in many studies the computational domain is limited to the wellbore zone or to the reservoir domain whereas the overpressure induced by the CO2 injection may extend several tens or hundreds kilometers. An artificial extension of the computational domain is often assumed to move away the boundary conditions but its simplified modeling cannot simulate complex fluid circulation at the basin scale. Several reasons justify this domain reduction. The geochemical impact is expected to be located in the storage area; the amount of data is mostly insufficient; the computational capacity of the current computers and their robustness to manage the heterogeneities of time scale and properties limit the reactive transport models. The current fluid flow simulators model the boundary condition with constant pressure or zero flux or analytical aquifer. For each face of the computational domain, a boundary condition, is imposed constant over time. However, when the induced overpressure extends beyond the limits of the studied domain, this type of boundary conditions does not satisfactorily model the transfer with the outer zone that impacts the simulation results of the reduced model.

4.1 Coupling Method Description

A coupling of two models based on two domains fitting into each other is proposed to overcome this issue (Fig. 16). One model (“model A”) only simulates fluid flow and CO2 dissolution in a domain that extends to the basin scale, while the other (“model B”), applied to a smaller, included domain that corresponds to the CO2 storage complex defined in Section 3, deals with the fluid-rock interaction.

thumbnail Figure 16

Method of coupling between a no reactive transport model applied to a domain at the basin scale (model-A) and a reactive transport model applied to a CO2 storage complex (model-B).

The coupling method consists in a one-way coupling between the model A and the model B. There is a transfer of information (pressure variation at the domain B limits) from the model A to the model B but there is no feedback from B to A (porosity and permeability variation of the model B due to fluid-rock interactions is not considered in the model A).

The workflow of the coupling process is (Fig. 16):

  • the no reactive transport model A is simulated over the whole studied period;

  • the reactive transport model B is simulated over the same time period with a pressure evolution over time as dynamic boundary conditions. This pressure is computed by the A simulation at dates defined in the data set of the A simulation.

4.2 Application to the France Nord Project: Description of the Models A and B

For the France Nord project, the domain (model) A corresponds to the entire domain described in Section 2 and the domain (model) B to the CO2 storage complex assessed in Section 3 (Fig. 17). However, the geochemical modeling is limited to the Chaunoy formation because it is the only reservoir layer for which enough chemical and mineralogical data are available. Therefore, the study of the geochemical effects focuses on reservoir sandstone formation.

thumbnail Figure 17

Chaunoy sandstone: initial pressure map. Dead cells (blue color) have no porosity and permeability values and are not included in the model.

The no reactive transport model applied to the domain A is the one detailed in Section 2.2 with a CO2 dissolution modeling. The reactive transport model applied to the domain B is the iterative coupling of the no reactive transport model A limited to the domain B with the geochemical model presented in Section 3.3. The permeability evolution due to porosity variation is not considered because the simulation results do not show a significant porosity variation.

The geochemical effects of CO2 injection are studied for only one well to simplify the local grid refinement in the injection area. It is located in the middle of the four wells defined in Section 3 (Fig. 17).

As explained in point 2 of the coupling workflow (Fig. 16), the user has to define a frequency of data transfer between the model A and the model B. For this case study, several runs were carried out to compare the injection rates and the gas saturation maps according to the frequency of the pressure outputs of the model A. The optimum frequency resulted to be annual, during the injection and the pressure stabilization period, then 10 years up to 150 years of storage, finally 50 years from 150 years to 1,000 years of storage.

4.3 Results

A comparison of the model A, the model B with static boundary conditions (named static BC) and the model B with dynamic boundary conditions (named dynamic BC), was carried out at 50 years to analyze the impact of each boundary condition on the pressure field (Fig. 18). The model with static BC is not coupled with the model A since the pressure at its boundaries is constant over time and fixed at the initial pressure. In the Figure 18a, the CO2 injection causes an overpressure which spreads out more than one hundred kilometers in the North-South direction therefore, beyond the boundaries of the model B domain. The pressure field of the model with dynamic BC is similar to the one of the model A. In the case of the model with static BC, an underestimation of the pressure is noticed. It is caused by an overestimation of the fluid flux at the limits of the domain B because the hydrostatic pressure as static BC is lower than the evolving pressure induced by the CO2 injection. With pressure as injection criterion, this leads to an overestimation of the injected CO2 mass of 10% compared to the dynamic case (Fig. 19) but not to a significant difference (Fig. 20) regarding the different trapping types (structural and stratigraphic, dissolution and mineral).

thumbnail Figure 18

Simulated pressure maps at 50 years: a) pressure map at the basin scale simulated by the model A, b) pressure map simulated by the model B with dynamic boundary conditions, c) pressure map simulated by the model B with static boundary conditions.

thumbnail Figure 19

Injected CO2 mass over time simulated by the model B: with static or dynamic boundary conditions.

thumbnail Figure 20

CO2 mass per trapping mechanism over time: static (SBC) or dynamic (DBC) boundary conditions.

In order to more precisely analyze the geochemical effects of the CO2 injection, several included local grid refinements were carried out: cell size from 3 km, 1 km, 300 m and to 100 m (Fig. 21).

thumbnail Figure 21

Refined mesh of the model B (storage complex).

Figures 22 and 23 show maps of gas saturation, dissolved CO2 and precipitated carbonates at 50 and 1,050 years. During the injection period, the CO2 advection caused by the induced overpressure is the main mechanism of transfer, so CO2 overall migrates in all the directions even if there is a preferential direction towards south-west, i.e., towards the Étampes fault because of the ascending slopes. The fault is not reached by the CO2 plume due to the structural trap between the fault and the injection well. With time, gravity becomes preponderant which explains the southwards CO2 (supercritical and dissolved) migration. Regarding the carbonation of the injected CO2, both the amount of carbonated CO2 and the carbonated mineral type evolve with time. At the beginning of the storage, CO2 is mineralized mainly as calcite and as dolomite at a negligible amount. With time, albite dissolves and allows dawsonite to precipitate to the detriment of calcite and dolomite that partly dissolve. At 1,050 years, dawsonite is the major carbonate. The slow kinetics of dawsonite compared to calcite causes a characteristic ring of calcite that surrounds the plume of dawsonite. The comparison of the CO2 plume (supercritical, dissolved, carbonates) location shows a difference between the plume of the fluid phase and the mineral phase one. Because of precipitation kinetics, the zone of precipitated carbonates corresponds to the dissolved CO2 plume 100 years ago. Regarding the porosity, no significant change is visible at 1,050 years because the rock volume lost by the dissolution of albite is compensated by the volume occupied by the dawsonite and quartz precipitation. In the dawsonite region, the porosity decreases of about 0.006 while it increases of 0.001 in the ring of calcite.

thumbnail Figure 22

Results of the model B coupled with the model A – Gas Saturation (SG) and CO2 molar fraction in brine (W-CO2s) at 50 years and 1,050 years. There is a visible effect of the coarse mesh where is based the refinement because the refinement uses a linear interpolation.

thumbnail Figure 23

Results of the model B coupled with the model A – mineral volume fraction of calcite (DPHIM-calcite) and of dawsonite (DPHIM-dawsonite) at 50 years and 1,050 years. There is a visible effect of the coarse mesh where is based the refinement because the refinement uses a linear interpolation.

In order to study the CO2 plume stabilization, the CO2 mass balance per trapping type is computed over time based on the simulation results. Figure 24 shows that during the 100 years of storage, CO2 is mainly supercritical, only 20% is dissolved in brine. From 100 years, the carbonates precipitate. At 1,050 years, about 78% of CO2 is carbonated mainly as dawsonite, 18% is dissolved and only 4% is supercritical. So, these results show a CO2 plume stabilization.

thumbnail Figure 24

CO2 mass distribution per trapping mechanism over time.

5 Discussion

The method, which consists in coupling two models (no reactive transport model A and reactive transport model B) based on domains fitting into each other was applied to test the injection of 200 Mt of CO2 in the Triassic formations of the South Keuper area. The simulation results have highlighted that the dynamic boundary conditions of the reactive transport model improve the fluid flow simulation and consequently the pressure field and the CO2 plume migration.

This application shows it would be relevant to use this approach in CO2 storage context. Indeed, for the performance assessment of a storage site, it is essential to estimate correctly its injectivity and its capacity for storing. In the same way, as part of a risk analysis, reducing the uncertainty of the pressure evolution in a CO2 storage site allows to decrease the uncertainty regarding the CO2 plume migration and the risk of fracturing the cap rock and of activating faults. Nevertheless, this proposed method needs to be improved because it is only one-way. Indeed, the impact of the porosity and permeability change at the storage-complex scale (model B) is not taken into account in the no reactive transport model (model A). For the France Nord project, the major fluid rock interactions occur after the injection period. At that time, the overpressure is reduced (lower than 2 MPa) and the porosity change is negligible. As a result, the one-way coupling is sufficient. But in case of quick and/or significant geochemical impact on porosity and permeability, the retroaction modeling would be necessary.

Regarding the CO2 mass balance per trapping type, the reactive transport simulation estimates a CO2 carbonation rate of about 78% at 1,050 years that seems to be excessive compared to the simulation study of other storage sites (Audigane et al., 2007; IPCC, 2005).

One of the possible causes of this excessive amount is the mesh refinement level. A sensitivity study using a coarse mesh (cell size of 3 km) and a fine mesh (cell size of 3 km) was carried out to estimate in a qualitative way the impact of the refinement on the CO2 behavior (Fig. 25). The numerical diffusion produced by the space discretization of the studied domain causes an overestimation of the CO2 dissolution. The direct consequence is a CO2 mineralization, which increases with the cell size. A comparison between Figures 24 and 26 show that at 1,050 years, (1) the rate of mineralized CO2 decreases from 88% to 78% thanks to the refinement, (2) there is still 4% of CO2 gas with the fine mesh whereas there is no gas anymore since 800 years with the coarse mesh.

thumbnail Figure 25

Volume fraction of dawsonite (DPHIM-dawsonite) at 1,050 years; a) coarse mesh (cell size of 3 km) and b) fine mesh (cell size of 100 m).

thumbnail Figure 26

CO2 mass distribution per trapping mechanism over time for a coarse mesh (3 km cell size).

However, this does not totally explain the excessive amount of carbonated CO2. The other main source of uncertainty concerns the geochemical model in particular, the secondary minerals and the precipitation kinetics (equilibrium constant, activation energy and the mineral surface where the geochemical reactions occur). An additional work dealing with both the kinetics data base and the textural models would be necessary in order to reduce the uncertainty of the injected CO2 mineralization.

Conclusion

The results presented were obtained in the framework of a preliminary desk study carried out before a possible CCS demonstration project in Upper Triassic formations of the Paris Basin. The objectives were to characterize a storage complex able to confine 200 Mt CO2 injected during 50 years, to simulate over time the CO2 migration and the induced pressure field, to optimize the number of wells, and to analyze the geochemical behavior of the porous media on the long term (1,000 years).

From previous work an appropriate area called “South Keuper sector” has been selected ca. 70 km south of Paris. It offers a potential reservoir in the Chaunoy sandstones, the thickness of which, ca. 60 m, was partly controlled by faults (the Pithiviers graben), an immediate caprock of at least local extension, the Rhaetian dolomitic Châlain red marls, and a regional 200 m thick confining formation, the Lower Liassic marine marls.

Under the major hypothesis that the reservoir formation behaves as a connected aquifer at the sector scale, the simulation of the base scenario showed that a ≥0.1 MPa overpressure area of 16,415 km2 concerns the Chaunoy formation at the end of the injection period (50 years), and spreads out to 51,869 km2 after 1,000 years. Overpressure peaks at 1.14 and 0.25 MPa, respectively, for the two stages. The large extension of pressure anomaly means that some brine migration at the outcrops could occur. Nevertheless the plume extension remains fairly small, respectively 316 km2 and 524 km2 at 50 and 1,000 years.

We simulated the fluid-rock interactions in a domain limited to the storage complex. For this purpose we developed a “one-way coupling” method between two numerical models. The first model simulates fluid flow and CO2 dissolution in a domain that covers the western part of Paris Basin. The second one deals with the fluid-rock interactions in a smaller, included domain that corresponds to the CO2 storage complex. The appropriate boundary conditions of the second model are derived from the results of the first model.

The results of the geochemical simulation showed that the injected CO2 tends to stabilize thanks to CO2 carbonation in calcite and dawsonite (dawsonite was considered, but a full solid solution model for clay minerals would be better adapted in a future work). The area affected by carbonation reaches an extension that depends on reaction kinetics. As dawsonite precipitates at a slower rate than calcite, a ring of calcite surrounding ring of dawsonite was depicted in the situation at 1,050 years. These fluid rock interactions did not cause any significant porosity change.

At the end of the simulation a CO2 carbonation rate of about 78% was calculated, that seemed somewhat high compared to simulations carried out on other storage sites. A sensitivity study discussed two possible reasons for such a difference: (1) the mesh refinement as source of excessive dissolution; (2) the geochemical model, in particular the too simple secondary minerals represented, and the precipitation kinetics considered.

The models built for this work are an essential contribution to the delimitation of the storage complex in the permeable formations. However, the geochemical and geomechanical impact on the cap rock remains to study to precise the vertical delimitation of the complex.

Acknowledgments

This work has been done within the France Nord project funded by ADEME (French Environment and Energy Management Agency) and industrial partners (TOTAL, ENGIE, EDF, Lafarge, Air Liquide, and Vallourec). The authors thank all the industrial partners in particular TOTAL and ENGIE for providing data. The authors are also grateful to N. Thybaud from ADEME, D. Copin from TOTAL and to all the partners – research institutes (BRGM, INERIS, and EFER) and industrial partners – for permission to publish these results. Finally, the authors thank all the technical contributors that have enabled to perform this work.

References

  • Audigane P., Gaus I., Czernichowski-Lauriol I., Pruess K., Xu T. (2007) Two-dimensional reactive transport modeling of CO2 injection in a saline aquifer at the Sleipner site, Am. J. Sci. 307, 2007, 974–1008. [CrossRef] [Google Scholar]
  • Bader A.G., Thibeau S., Vincké O., Delprat Jannaud F., Saysset S., Joffre G.H., Giger F.M., David M., Gimenez M., Dieulin A., Copin D. (2014) CO2 storage capacity evaluation in deep saline aquifers for an industrial pilot selection. Methodology and results of the France Nord project, Energy Procedia 63, 2779–2788. [CrossRef] [Google Scholar]
  • Baker J.C., Bai G.P., Hamilton P.J., Golding S.D., Keene J.B. (1995) Continental-scale magmatic carbon dioxide seepage recorded by dawsonite in the Bowen-Gunnedah-Sydney basin system, Eastern Australia, Int. J. Sediment Res. A65, 3, 522–530. [Google Scholar]
  • Bazin B., Brosse É, Sommer F. (1997) Chemistry of oil-field brines in relation to diagenesis of reservoirs. 1. Use of mineral stability fields to reconstruct in situ water composition. Example of the Mahakam basin, Mar. Pet. Geol. 14, 5, 481–495. [CrossRef] [Google Scholar]
  • Bénézeth P., Palmer D.A., Anovitz L.M., Horita J. (2007) Dawsonite synthesis and reevaluation of its thermodynamic properties from solubility measurements: Implications for mineral trapping of CO2, Geochim. Cosmochim. Acta 71, 4438–4455. [CrossRef] [Google Scholar]
  • Bennehard M. (1970) Contribution à l’étude minéralogique et sédimentologique des terrains triasiques dans quelques sondages du bassin de Paris, PhD Univ, Toulouse, p. 282. [Google Scholar]
  • Bouc O., Fabriol H., Brosse É., Kalaydjian F., Farret R., Gombert P.H., Berest P., Lagneau V., Pereira J.-M., Fen-Chong T. (2011) Lignes de conduite pour la sécurité d’un site de stockage géologique de CO2, BRGM/RP-60369-FR 154, 3 annexes. [EDP Sciences] [Google Scholar]
  • Bourquin S., Robin C., Guillocheau F., Gaulier J.-M. (2002) Three-dimensional analysis of the Keuper of the Paris Basin: Discrimination between tectonics, eustasy and sediment supply in the stratigraphic record, Mar. Pet. Geol. 19, 469–498. [Google Scholar]
  • BRGM (2008) Évaluation du potentiel géothermique des réservoirs clastiques du Trias du Bassin de Paris, Rapport final RP-56463-FR. [Google Scholar]
  • Brooks R.H., Corey A.T. (1966) Properties of porous media affecting fluid flow, J. Irrig. Drain. Eng., 92, 61–90. [Google Scholar]
  • Brosse E., Fabriol H., Fleury M., Grataloup S., Lombard J.M. (2010a) CO2 storage in the struggle against climate change, Oil Gas Sci. Technol. – Rev. IFP 65, 3, 369–373. [CrossRef] [EDP Sciences] [Google Scholar]
  • Brosse É., Badinier G., Blanchard F., Caspard E., Collin P.Y., Delmas J., Vidal-Gilbert S. (2010b) Selection and characterization of geological sites able to host a pilot-scale CO2 storage in the Paris Basin (GéoCarbone-PICOREF), Oil Gas Sci. Technol. – Rev. IFP 65, 3, 375–403. [Google Scholar]
  • Comité des Techniciens ELF (1991) Monographies des principaux champs pétroliers de France. Bulletin des centres de recherches Exploration-production Elf-Aquitaine. Mémoire 14. 2-901026-354. [Google Scholar]
  • Coudrain-Ribstein A., Gouze P., de Marsily G. (1998) Temperature-carbone dioxide partial pressure trends in confined aquifers, Chem. Geol. 145, 73–89. [CrossRef] [Google Scholar]
  • Delmas J., Brosse É., Houel P. (2010) Petrophysical properties of the middle Jurassic carbonates in the PICOREF Sector (South Champagne, Paris Basin, France), Oil Gas Sci. Technol. – Rev. IFP 65, 3, 405–434. [Google Scholar]
  • Delmas J., Houel P., Vially R. (2002) Paris Basin – petroleum potential, IFP Regional Report n° 59994. [Google Scholar]
  • De Lucia M., Kempka T., Kühn M. (2014) A coupling alternative to reactive transport simulations for long-term prediction of chemical reactions in heterogeneous CO2 storage systems, Geosci. Model Dev. Discuss. 7, 6217–6261. [CrossRef] [Google Scholar]
  • Eschard R. et al. (1998) Combining sequence stratigraphy, geostatistical simulations, and production data for modeling a fluvial reservoir in the Chaunoy field (Triassic, France), AAPG Bull. 32, 4, 545–568. [Google Scholar]
  • Fornel A., Roggero F., Vincké O. (2011) Simulations d’injection de CO2 dans l’aquifère salin du Trias Ouest – Projet France Nord, IFPEN, Report No. 61765. [Google Scholar]
  • Fornel A., Delaplace P., Vincké O. (2012) Projet France Nord - Sensibilité des paramètres du modèle sur la capacité de stockage de CO2 en 50 ans dans la zone du Keuper Sud, IFPEN, Report No. 62575. [Google Scholar]
  • Gao Y., Liu L., Hu W. (2009) Petrology and isotopic geochemistry of dawsonite-bearing sandstones in Hailaer basin, northeastern China, Appl. Geochem. 24, 1724–1738. [CrossRef] [Google Scholar]
  • Gonçalvès J., Violette S., Guillocheau F., Robin C., Pagel M., Bruel D., Ledoux E. (2004) Contribution of a three-dimensional regional scale basin model to the study of the past fluid flow evolution and the present hydrology of the Paris Basin, France, Basin Res. 16, 4, 569–586. [CrossRef] [Google Scholar]
  • Guillocheau F., Robin C., Allemand P., Bourquin S., Brault N., Dromart G., Friedenberg R., Garcia J.-P., Gaulier J.-M., Gaumet F., Grosdoy B., Hanot F., Le Strat P., Mettraux M., Nalpas T., Prijac C., Rigoltet C., Serrano O., Grandjean G. (2000) Meso-Cenozoic geodynamic evolution of the Paris Basin: 3D stratigraphic constraints, Geodin. Acta 13, 4, 189–245. [Google Scholar]
  • Hellevang H., Declercq J., Aagaard P. (2011) Why is dawsonite absent in CO2 charged reservoirs? Oil Gas Sci. Technol. – Rev. IFP 66, 1, 119–135. [CrossRef] [EDP Sciences] [Google Scholar]
  • Houel P., Euzen T. (1999) Atlas Régional du Trias de l’Ouest du Bassin Parisien, (modélisation des aquifères), Rapport IFP No. 53171 [Google Scholar]
  • Intergovernmental Panel on Climate Change (2005) IPCC special report on carbon dioxide capture and storage, Cambridge University Press, New York, NY, USA. [Google Scholar]
  • Johnson J.W., Nitao J.J., Knauss K.G. (2004) Reactive transport modeling of CO2 storage in saline aquifers to elucidate fundamental processes, trapping mechanisms and sequestration partitioning, in: S.J. Bains, R.H. Worden (eds), Geological storage of carbon dioxide, Geological Society Special Publications, London, pp. 107–128. [Google Scholar]
  • Kaszuba J.P., Viswanathan H.S., Carey J.W. (2011) Relative stability and significance of dawsonite and aluminum minerals in geologic carbon sequestration, Geophys. Res. Lett. 38, L08404. [CrossRef] [Google Scholar]
  • Kharaka Y.K., Hanor J.S. (2003) Deep fluids in the continents: I. Sedimentary basins, in: Drever J.I. (ed.), Holland H.D., Turekian K.K. (Executive Editors), Treatise on Geochemistry, Vol. 5, Elsever, pp. 499–540, ISBN 0-08-043751-6. [Google Scholar]
  • Lasaga A.C. (1984) Chemical kinetics of water-rock interactions, J. Geophys. Res. B 89, B6, 4009–4025. [CrossRef] [Google Scholar]
  • Lasaga A.C. (1995) Fundamental approaches to describing mineral dissolution and precipitation rates, in: White A.F., Brantley S.L. (eds), Reviews in mineralogy volume 31: Chemical weathering rates of silicate minerals, Washington, DC, Mineralogical Society of America, pp. 23–86. [Google Scholar]
  • Lasaga A.C. (1998) Kinetic theory in Earth Sciences, Princeton University Press, Princeton. [Google Scholar]
  • Lecomte J.-C., Houel P., Daniel J-M., Vincké O. (2010) Modélisation géologique du Trias de l’Ouest du Bassin de Paris. Sélection préliminaire de zones d’injection et de stockage de CO2 pour le projet France-Nord, IFPEN, Report No. 61670. [Google Scholar]
  • Le Gallo Y., Trenty L., Michel A., Vidal-Gilbert S., Parra T., Jeannin L. (2006) Long-term flow simulation of CO2 storage in saline aquifers, in Proceedings of the 8th International Conference on Greenhouse Gas Control Technologies, 19-23 June, Trondheim, Norway. [Google Scholar]
  • Lohrenz J., Bray B., Clark C. (1964) Calculating viscosities of reservoir fluids from their compositions, J. Pet. Technol. 16, 1171–1176. [CrossRef] [Google Scholar]
  • Moutte J. (2009) Arxim, a library for thermodynamic modeling of fluid-rock systems, 55, (http://www.emse.fr/~moutte/arxim/). [Google Scholar]
  • Moutte J., Michel A., Battaia G., Parra T., Garcia D., Wolf S. (2010) Arxim, a library for thermodynamic modeling of reactive heterogeneous systems, with applications to the simulation of fluid-rock systems, 21st Congress of IUPAC, Conference on Chemical Thermodynamics, Tsukuba, Japan. [Google Scholar]
  • Newman G.H. (1973) Pore-volume compressibility of consolidated, friable, and unconsolidated reservoir rocks under hydrostatic loading, J. Pet. Technol. 25, 2, 129–134. [CrossRef] [Google Scholar]
  • Okuyama Y., Todaka N., Sasaki M., Ajima S., Akasaka C. (2013) Reactive transport simulation study of geochemical CO2 trapping on the Tokyo Bay model – With focus on the behavior of dawsonite, Appl. Geochem. 30, 57–66. [CrossRef] [Google Scholar]
  • Palandri J.L., Kharaka Y.K. (2004) A compilation of rate parameters of water-mineral interaction kinetics for application to geochemical modeling, U.S.G.S Open File Report 2004-1068, 64. [EDP Sciences] [Google Scholar]
  • Peng D.Y., Robinson D.B. (1976) A new two-constant equation of state, Ind. Eng. Chem. Fundam. 15, 59–64. [CrossRef] [Google Scholar]
  • Perrette L. (2011) Transposition of the CCS Directive in France, IEG Meeting, 14 November, Brussels, 22, Presentation. [Google Scholar]
  • Perrodon A., Zabek J. (1990) Paris Basin. Interior cratonic basins, AAPG Memoir n° 51, pp. 633–679. [Google Scholar]
  • Reid R.C., Prausnitz J.M., Poling B.E. (1987) The Properties of Gases and Liquids, 4th edn., McGraw-Hill, New York. [Google Scholar]
  • Schmidt E. (1969) Properties of water and steam in SI units, Springer Verlag, Berlin. [Google Scholar]
  • Serra O. (1985) Diagraphies différées. Bases de l’interprétation. Tome 2. Interprétation des données diagraphiques. Bulletin des centres de recherches exploration-production Elf Aquitaine. Memoire 7, ISSN 0181-0901. [Google Scholar]
  • Spötl C., Wright V.P. (1992) Groundwater dolocretes from the Upper Triassic of the Paris Basin, France: a case study of an arid, continental diagenetic facies, Sedimentology 39, 1119–1136. [CrossRef] [Google Scholar]
  • Soreide I., Whitson C.H. (1992) Peng-Robinson predictions for hydrocarbons, CO2, N2, and H2S with pure water and NaCl brine, Fluid Phase Equilib. 77, 217–240. [CrossRef] [Google Scholar]
  • Thomas L.K., Katz D.L., Tek M.R. (1968) Threshold pressure phenomena in porous media, SPE J 8, 174–184. [CrossRef] [Google Scholar]
  • Wertz F., Gherardi F., Blanc P., Bader A.-G., Fabbri A. (2012) Modelling CO2-driven cement alteration at well-caprock interface, Proceedings, Tough Symposium 2012 Lawrence Berkeley National Laboratory, 17-19 September, Berkeley, California. [Google Scholar]

All Tables

Table 1

Mean value of the properties for each geological formation modeled for this study.

Table 2

CO2 thermodynamic properties (Reid et al., 1987).

Table 3

Water analysis provided for test N.2 at well Beaumont N.101 (Triassic layers ca. 2,240 m logging depth – close to true vertical depth).

Table 4

Composition and kinetic parameters considered for representing the reservoir in the geochemical simulations.

Table 5

Impact of the injection on exploited oil and gas fields.

All Figures

thumbnail Figure 1

Stratigraphic column of the South Keuper sector area.

In the text
thumbnail Figure 2

Structural map at the top of the Triassic formations. The potential CO2 storage sites should be located in the areas represented by either the black rectangle (South Keuper) or the yellow one (North Keuper). The orange meridian line is the eastern limit of the geological model.

In the text
thumbnail Figure 3

Hydraulic layers. The South Keuper sector area is located in the yellow rectangle.

In the text
thumbnail Figure 4

“Quick-Look” layering of the Norian Chaunoy sandstones formation.

In the text
thumbnail Figure 5

Mesh of the geological model and schematic view of the boundary conditions, for the selected area: imposed pressure at the outcrops (0.1 MPa at 200 m subsea, in green, and 0.1 MPa at 0 m, in purple). The mesh represents the top of the Triassic formations. It is limited at the eastern boundary by the meridian line defined in Figure 2.

In the text
thumbnail Figure 6

Structural map showing faults and the Triassic top layer: a) faults (black and pink lines) of the Paris Basin; b) faults of the geological model illustrated in Figure 5.

In the text
thumbnail Figure 7

Triassic outcrops (purple area) at the west of the Paris Basin. Overlay of cultural data, surface topography (radar image) and isobaths (m) at the top of the Triassic series.

In the text
thumbnail Figure 8

Structural grid with porosity values and location of the injection well defined in Section 3.1.

In the text
thumbnail Figure 9

Water and gas a) relative permeability curves and b) capillary pressure (0.1 MPa) curve.

In the text
thumbnail Figure 10

a) Location of the control wells chosen in the Paris Basin in order to b) compare measured pressure values, P gauge, and simulated ones, P sim (right).

In the text
thumbnail Figure 11

Mass of injected CO2.

In the text
thumbnail Figure 12

Maximal overpressure (MPa) in the Chaunoy formation over time (year).

In the text
thumbnail Figure 13

Overpressure (MPa) above 0.1 MPa at the a) end of the injection and b) 1,000 years after – Chaunoy layer.

In the text
thumbnail Figure 14

CO2 plume (gas saturation) at the a) end of the injection and b) 1,000 years after – Chaunoy layer. The limits of the storage complex are represented by the red rectangle.

In the text
thumbnail Figure 15

Sensitivities (red) and induced errors (green) of uncertain parameters on the mass of injectable CO2. (1) Maximum overpressure at injection point (uncertainty of 1% of the initial pressure), (2) porosity (uncertainty of 50% of the initial values, (3) horizontal permeability (uncertainty of 80% of the initial values), (4) rock compressibility (uncertainty of 50% of the initial value), (5) faults sealing (multiplying transversal transmissibility factor equal to zero or 1), (6) CO2 dissolution (modeled or not), (7) irreducible water saturation (uncertainty of 50% of the initial value) with dissolution modeling, (8) vertical permeability anisotropy ratio (0.1 or 0.0001), (9) irreducible water saturation (uncertainty of 50% of the initial value) with no dissolution modeling.

In the text
thumbnail Figure 16

Method of coupling between a no reactive transport model applied to a domain at the basin scale (model-A) and a reactive transport model applied to a CO2 storage complex (model-B).

In the text
thumbnail Figure 17

Chaunoy sandstone: initial pressure map. Dead cells (blue color) have no porosity and permeability values and are not included in the model.

In the text
thumbnail Figure 18

Simulated pressure maps at 50 years: a) pressure map at the basin scale simulated by the model A, b) pressure map simulated by the model B with dynamic boundary conditions, c) pressure map simulated by the model B with static boundary conditions.

In the text
thumbnail Figure 19

Injected CO2 mass over time simulated by the model B: with static or dynamic boundary conditions.

In the text
thumbnail Figure 20

CO2 mass per trapping mechanism over time: static (SBC) or dynamic (DBC) boundary conditions.

In the text
thumbnail Figure 21

Refined mesh of the model B (storage complex).

In the text
thumbnail Figure 22

Results of the model B coupled with the model A – Gas Saturation (SG) and CO2 molar fraction in brine (W-CO2s) at 50 years and 1,050 years. There is a visible effect of the coarse mesh where is based the refinement because the refinement uses a linear interpolation.

In the text
thumbnail Figure 23

Results of the model B coupled with the model A – mineral volume fraction of calcite (DPHIM-calcite) and of dawsonite (DPHIM-dawsonite) at 50 years and 1,050 years. There is a visible effect of the coarse mesh where is based the refinement because the refinement uses a linear interpolation.

In the text
thumbnail Figure 24

CO2 mass distribution per trapping mechanism over time.

In the text
thumbnail Figure 25

Volume fraction of dawsonite (DPHIM-dawsonite) at 1,050 years; a) coarse mesh (cell size of 3 km) and b) fine mesh (cell size of 100 m).

In the text
thumbnail Figure 26

CO2 mass distribution per trapping mechanism over time for a coarse mesh (3 km cell size).

In the text

Les statistiques affichées correspondent au cumul d'une part des vues des résumés de l'article et d'autre part des vues et téléchargements de l'article plein-texte (PDF, Full-HTML, ePub... selon les formats disponibles) sur la platefome Vision4Press.

Les statistiques sont disponibles avec un délai de 48 à 96 heures et sont mises à jour quotidiennement en semaine.

Le chargement des statistiques peut être long.