Determination of the effect of geological reservoir variability on carbon dioxide storage using numerical experiments

Résumé — Détermination de la variabilité des réservoirs géologiques sur le stockage du CO 2 par la méthodologie des plans d’expériences — Dans le contexte de l’étude du stockage géologique du dioxyde de carbone dans les réservoirs sédimentaires, les simulations tentent de reproduire le plus fidèlement possible les écoulements fluides et gaz, ainsi que les réactions chimiques entre les minéraux (calcite et dolomite) et le CO 2 injecté [André et al . (2007) Energy Convers. Manage . 48, 1782-1797; Gunter et al . (1999) Appl. Geochem . 4, 1-11]. Cependant, les réservoirs étudiés sont généralement mal connus. Les observations sont peu nombreuses et par suite les variogrammes des propriétés pétrophysiques et minéralogiques sont approchés. L’article à incertitudes prévisions stockage . en incertitudes deux paramètres opérationnels : quantité dioxyde de carbone stockée variation moyenne porosité. retenues sont variabilité du tirage Abstract — Determination of the Effect of Geological Reservoir Variability on Carbon Dioxide Storage Using Numerical Experiments — The simulations of carbon dioxide storage in sedimentary reservoirs model the fluid and gas flow and the chemical reactions which occur between the minerals and dolomite) and the injected CO 2 et al . (2007) Energy Convers. 1782-1797; . data, focus two operational parameters: the quantity of the stored carbon dioxide and the mean variation of the porosity. Two sources of uncertainties are examined: the draw dispersion and the approximation on the variogram parameters. To study the influence of the draw dispersion, variogram parameters are kept fixed and different simulations are run; the associated variance on the operational parameters then has the meaning of a repeatability error. In the second a sensibility analysis is carried out to study the influence of variogram parameters variations (sill, range, nugget effect) on the CO 2 storage. The chosen methodology is the designs of experiments.


INTRODUCTION
With increasing international concern over global climate change probably due to increased CO 2 emissions, CO 2 capture and storage is emerging as a potentially important way of reducing CO 2 concentration in the atmosphere. Geological formations must have suitable petrophysical properties and well-adapted chemical composition in order to be good candidate media for CO 2 sequestration. The understandings of the geology of geological reservoirs and of the evolution of their petrophysical properties are then crucial. Especially a good appreciation of the different parameters controlling the reactions of precipitation and dissolution of minerals in porous media during the injection of CO 2 is necessary. Indeed, when an acid gas like CO 2 is injected in permeable rocks (saline aquifer or depleted hydrocarbon reservoir), this gas dissolves in the water present in the aquifer, which becomes acid, aggressive and dissolves some of the minerals (Lindeberg et al., 2002). This initiates a sequence of chemical reactions that can last for thousands of years (André et al., 2007).
In this context, numerical simulations through reactive transport software are available to model the complex phenomena that occur during the CO 2 injection. Simulations of the CO 2 storage in sedimentary reservoirs predict the behaviour of these reservoirs through the modeling of the fluid transfers and of the chemical reactions with the minerals (Brosse et al., 2005;Johnson et al., 2004); they also permit to evaluate possible mechanical consequences of this storage and the associated risks.
However, theses numerical simulations assume a reasonably complete understanding of the reservoir before the injection. But these reservoirs are often heterogeneous; their geological properties vary with space and evolve with the storage itself. Therefore, a realistic modeling of all geochemical reactions is not sufficient to clearly describe the evolution that occurs in the reservoir. There is a need for quantifying the effects of spatial heterogeneity and geological properties uncertainties on the result of the CO 2 sequestration (Chilès and Delfiner, 1999). Indeed, the spatial initial variability of the different reservoirs characteristics (porosity and permeability of the rock, proportions of various minerals, etc.) has a major effect on the behavior of the system, because of the creation of preferential fluid flows that delay or increase the advancement of the reaction front (De Lucia, 2008).
The aim of the paper is to quantify the impact of the initial variability of the geological reservoir on the underground storage of carbon dioxide (Hovorka et al., 2004). The approach is similar to (De Lucia et al., 2011) study but with a more complex geochemical system to be closer to reality (Michard and Bastide, 1988). In this study, a carbonated reservoir from the Dogger of Paris Basin is considered

Abstract -Determination of the Effect of Geological Reservoir Variability on Carbon Dioxide
Storage Using Numerical Experiments -The simulations of carbon dioxide storage in sedimentary reservoirs model the fluid and gas flow and the chemical reactions which occur between the minerals (calcite and dolomite) and the injected CO 2 [André et al. (2007) Energy Convers. Manage. 48, 1782-1797Gunter et al. (1999) Appl. Geochem. 4,[1][2][3][4][5][6][7][8][9][10][11]. However because of the lack of data, these reservoirs are always partially known and the fitted variograms of petrophysical and mineralogical quantities are approximate. The aim is to quantify the impact of uncertainties on reservoir characteristics on the storage predictions. We focus on two operational parameters: the quantity of the stored carbon dioxide and the mean variation of the porosity. Two sources of uncertainties are examined: the draw dispersion and the approximation on the variogram parameters. To study the influence of the draw dispersion, variogram parameters are kept fixed and different simulations are run; the associated variance on the operational parameters then has the meaning of a repeatability error. In the second case, a sensibility analysis is carried out to study the influence of variogram parameters variations (sill, range, nugget effect) on the CO 2 storage. The chosen methodology is the designs of experiments. The simulations are carried out using reactive transport software. The studied carbonated reservoir is built in reference to the Dogger formation of Paris Basin (France) [Diedro (2009) (Delmas, 2007;Brosse et al., 2007, Rojas et al., 1989. We focus on the Dogger formation because it is expected to store carbon dioxide in carbonated reservoir similar to the Dalle Nacrée, the white Oolith and the Comblanchien.
The CO 2 storage is examined through two operational parameters evolving with time: the quantity of carbon dioxide stored and the mean variation of the porosity in the reservoir.
Two sources of uncertainties are studied: the draw dispersion and the variations of the fitted variogram parameters. In the first case, simulations are run on different realizations of the same random field (same variogram parameters). The observed variance on the operational parameters has the meaning of a repeatability error. In the second case, we study the effects of mineralogical and petrophysical parameter variations. The considered variations are induced by perturbations of the variogram parameters (sill, range, nugget effect). The aim is then to measure the consequences of a possible "variographic fitting error". The methodology of Design of Experiments (DOE) is used to carry out the sensibility analysis. The principle consists in varying simultaneously the levels of several factors within each test (Helbert, 2006;Linder, 2005). The influence of some parameters and their interaction is detected using a minimal number of simulations while keeping and preserving the precision of the results.
The simulations are carried out using the Coores software of reactive transport (Diedro, 2009). Coores is developed in partnership between the laboratory of geochemistry of the École Nationale Supérieure des Mines de St-Étienne and with IFP Energies nouvelles. It couples a geochemistry module and a transport module. The transport module developed at IFP Energies nouvelles calculates the displacements of fluid (water and gas) following the generalized law of Darcy. The geochemistry module developed at the École des Mines computes the modifications due to the composition of the fluid and evaluates their consequences on the petrophysical properties (Michel and Trenty, 2005).
The basic equation for reactive transport modelling states the conservation of the chemical species in the aqueous phase. For a given volume V tot of the porous medium, let n aq l be the number of moles of species l. This number may vary because of the addition or subtraction of the chemical species by means of the pervading fluid through the rock (the velocity of the fluid is u and the porosity of the rock is Φ). It may also change because of the dissolution or precipitation of the minerals m contained in the solid portion of the rock. Let Φ m be the proportion of solid phase occupied by mineral m, the molar volume of which is V m . And let β l,m be the number of moles that mineral m may add to or subtract from the aqueous solution depending on its dissolution or precipitation. The conservation equation then reads: This equation has to be coupled with the laws ruling the kinetics of the mineral dissolution or precipitation, as a function of the thermodynamic disequilibrium expressed in the ratio Q m /K m ; K m is the equilibrium constant for mass action law and Q m is the activity product, both parameters corresponding to the chemical reaction written for the dissolution/precipitation of mineral m in water. The kinetic law is also a function of the kinetic constant k m and of the surface S m of the mineral (a factor depending on pH may also be added).
Appropriate chemical, thermodynamic and kinetic data for the different minerals must be specified to the computer code (Gaucher et al., 2006). The velocity law is dependent on Darcy's law and is linked to the pressure and porosity/permeability field (Aubertin et al., 1987). All the parameters, whatever chemical, thermodynamic, petrophysical or related to fluid dynamics could be subject to variation. As regards the influence on CO 2 storage, only the variations of the mole numbers of the minerals are examined and the variations of the other parameters are disregarded.
Section 1 introduces the modelling of the carbonated reservoir in reference to the Paris Basin and the simulated system. In Section 2, we present the results of a preliminary study carried out to quantify the impact of the draw variability for a fixed set of variogram parameters. The impact of variogram parameters on the CO 2 storage is considered in Sections 3 and 4: the methodology of experimental design is presented in Section 3, the mathematical model is then evaluated and the results are analysed in Section 4.
As an example, Figure 1 displays the output obtained by software Coores when CO 2 is injected in the reservoir: the residual gas is progressively dissolved in the water of the aquifer (Diedro, 2009;Castro et al., 1998a,b).

SELECTED RESERVOIRS AND REACTIVE
TRANSPORT SIMULATIONS

Modeling of the Carbonated Reservoir in Reference to the Paris Basin
The simulations should fit some similarity with real geological instances; thus, they were inspired by the levels from the Dogger of Paris Basin that are judged as a candidate for CO 2 storage. The selected levels are the "Oolithe Blanche", the "Comblanchien" (upper Bathonian) and the Dalle Nacrée (lower Callovian). The simulated model is based on the Dalle Nacrée, in order to study the influence of the clay variability (clay content about 17%; Kohler, 2006) and the Comblanchien and the Oolithe Blanche, containing less than one percent clay. Petrophysical information comes from the St. Martin de Bossenay 1 drill core (SMB1) in the Oolithe Blanche (between 1 436.35 m and 1 542 m; Delmas, 2007). The mineralogical information comes from the unique core (LCV) from Ravières. Other information comes from four drills in the Oolithe Blanche and three in the Comblanchien. The information on the nature and proportion of minerals has been obtained by combining petrographical information and DRX (X-ray diffraction) measurements on the rock samples taken from the drill cores. In some cases, use of electron microprobe has been made in order to know the exact chemical composition of the minerals. The more precise proportion of minerals has been obtained by using the data on whole rock geochemistry and by converting them to mineral modes, thanks to the information on the nature and composition of the present minerals.
The petrophysical and mineralogical variability (variogram, means and variance) of the fields also have been obtained from the geochemical analysis of verticals wholes rock. Then, multivariate geostatistics methods permit to stimulate correlated variables, such as the couple porosity, permeability and also the minerals. Finally a multiplicative coefficient was applied on the range and permitted the passage of the vertical field to horizontal field in two dimensions.
The main minerals observed are silicates (quartz, illite, chlorite, kaolinite) and carbonates (calcite and dolomite) (Azaroual et al., 1997). Field observations show that dolomite is not present uniformly at the decimeter to meter scale but forms islands or disjoint spots.
For the purpose of geochemical modeling, the silicates are considered as a single group noted TSM (total of silicate minerals) with quartz = 32.9% TSM, chlorite = 11% TSM, illite = 47.9% TSM, kaolinite = 8.2% TSM. The field of dolomite spots and the field of calcite are deduced from the closing relation, calculated by subtracting the other minerals. The vertical variogram is obtained using drilling information from Ravières (LCV) core. Since the horizontal variogram is impossible to obtain because of the lack of data, it is simply Diminution and dissolution of the CO 2 residual gas in the reservoir following the evolution of the aquifer from west to east. a) t = 500 years, b) t = 1 000 years, c) t = 1 500 years, d) t = 2 000 years.
derived from the vertical one, applying a proportional factor on the range; the factor is such that 10 m in the vertical field represent 300 m in the horizontal field, as suggested by several observations in the similar formations from Bahamas Dogger (Collin, 2007). The three groups of parameters give rise to eight factors: -the variogram of the logarithm of the dolomite spot field is characterized by two parameters: sill and range; -the silicate variability is characterized by three parameters: sill, range and nugget effect of Total Silicate Minerals (TSM) field; -petrophysic is characterized by three parameters: sill, range and nugget effect, for both porosity and permeability that are linked.

The Simulated System
The simulations were run until 2 000 years on a field of dimensions 1 000 m × 1 000 m, with a meshing 80 × 80 and with 12.5 m width. The considered reservoir is essentially composed of carbonates, mainly calcite, dolomite and silicate (illite, chlorite, kaolinite and quartz). The CO 2 was injected in a heterogeneous carbonated reservoir in residual saturation in the porous media at a pressure of 150 bar (Coudrain-Ribstein et al., 1998). A proportion of CO 2 is dissolved in the aquifer while the remaining stays in gaseous form and occupies 1/10 of pores volume. We consider the geological reservoir is a mobile aquifer moving at a velocity of 1 m/year (Aubertin et al., 1987;Palandri and Kharaka, 2004;Pruess et al., 2001). The contact between the aquifer and the gas modifies the equilibrium of the system and initiates the chemical reactions in the reservoir.

UNCERTAINTY OF THE INITIAL STATE
Variograms are needed for each parameter that can potentially influence the CO 2 storage. These quantities are petrophysical properties (porosity and permeability), silicates and dolomite spots within the rock. Three empirical variograms are thus calculated and the fitted model is the sum of a nugget component and a spherical component: ( 2) where: -Pep(h) = 1 if h ≠ 0 and 0 otherwise, λ 1 is the total sill, λ 2 is the nugget effect and λ 3 is the range.
Once the variograms parameters are fitted, a stochastic realization of the geological spatial field at the initial instant is obtained (Fig. 2).
This realization is then used as input of the reactive transport software. The result of the run is the evolution of the geological field with time. We particularly focus on the evolution of the quantity of carbon dioxide stored and the mean variation of the porosity in the reservoir.
To analyze the sensibility of the storage process to the geological field, two sources of uncertainties are examined. In the first one, all variograms parameters are fixed and different realizations of the initial field are used (Chilès and Delfiner, 1999). The impact on the storage of this variability is called the draw dispersion. In the second one, the parameters of the variogram are not fixed anymore. Different realizations of Observation of a stochastic realization: a) porosity geological field at time t = 0, b) variogram of the porosity field. the initial geological field for different values of range and sill are studied. Indeed, the parameters of the variogram (range, sill and nugget effect) are fitted on the available data. They can be imprecise because of the lack of data or measurement errors. Moreover, it is necessary to take account of the scaling between the grid and the measured values of the drilling. Therefore, the impact of variations of the variogram parameters around the reference values on carbon dioxide storage is examined.
Section 3 introduces the reference values of the parameters of interest calibrated on observations and also their field of variation.

VARIOGRAM CALIBRATION AND FIELD VARIATION OF PARAMETERS
The factors of which we want to quantify the impact on the output are parameters of the variograms that characterize the geological reservoir. The estimation of these parameters and their range of variation for silicates, petrophysical properties and for dolomite spots within the rock is introduced here.

Silicates
The Dalle Nacrée part of the field contains about 17% of clays mineral whereas the Oolith and the Comblanchien contain less than 1% (Kohler, 2006). We make the hypothesis that the variation of silicate is the same in the three groups (Dalle Nacrée, white Oolith and Comblanchien). The analysis shows that the concentrations of the silicate minerals are all proportional and are thus deducted from the Total Silicate Minerals (TSM) from the DRX (Kohler, 2006), referred as TSM. Their compositions have the same variability in those total silicate minerals (Fig. 3). The illite represents 32% of TSM, chlorite 11%, quartz 32.9% and kaolinite 8.2%. The vertical variogram fitted on the data of the white Oolith formation has the following parameters (see Eq. 2): the range λ 1 is equal to 0.16, the nugget effect λ 2 to 0.08 and the sill λ 3 to 12. The proposed domain of variation for λ 1 is [0.12; 0.2], for λ 2 [0.05; 0.11] and for λ 3 [6; 18]. Theses perturbation intervals are consistent with observations (Diedro, 2009;Collin, 2007).

Dolomite Spots Within the Rock
The geological and petrographical observations in the field show that the dolomite mineral is not present everywhere in the reservoir at the scale the rock is described. On the contrary, it can be considered that dolomite appears as spots at the macroscopic scale. We speak of the size of the dolomite spots within the rock to characterize this scattering. It is important to notice that we don't used here the usual Gaussian threshold modal but the presence of the dolomite spots is simulated using a threshold of the logarithm of the dolomite concentration field: so at fixed dolomite content, thanks to the concentrations variations, the greater the size of the spots, the larger the volume and the reactive surface, for a given size of mineral grains. Here to better fit the observation of the dolomite field, we choose these parameters of the variogram of the logarithmic transformation of the dolomite concentration: a range equal to 7 and a sill equal to 0.01. These values were obtained by the method of errors and testing, with a mean = 0.03 fitted on the data (Diedro, 2009). The variogram is defined by: where λ dol/1 is the sill and λ dol/2 is the range. The domain of variation of these quantities ([2, 12] for the range and [0.006, 0.014] for the sill) is chosen to allow small variations below and above the mean of the fitted values. The form of the dolomite field varies according to the range, to the sill and to the threshold. When the range is high, the threshold concentration is fixed to 0.019 so that it generates large spots of dolomite. In the other hand, when the range is low, the threshold is fixed to 0.033 and generates small spots of dolomite within the rock.
Four different fields are generated with the same random values. As illustrated in Figures 5 and 6, they represent the diversity of the size of the spots of dolomite within the rock.  Examples of porosity variograms. a) Low nugget effect (vario 1: small sill and small range; vario 2: large sill and small range; vario 3: small sill and high range; vario 4: high sill and high range). b) High nugget effect (vario 1: small sill and small range; vario 2: large sill and small range; vario 3: small sill and high range; vario 4: high sill and high range).

Petrophysics
We can observe on experimental data (Delmas, 2007) a link between porosity and permeability that can be modeled by a linear model of coregionalisation. The simple and the cross structural variogram define a linear model of the same basis structure. These two quantities can be jointly simulated. We use the following variogram of the transformed Gaussian with a linear model of coregionalisation. Thus, we only study the impact of the parameters of the variogram of porosity. The chosen variogram is the following: (4) where λ poro1 represents the total range, λ poro2 the nugget effect and λ poro3 the sill. The experimental domain is calibrated using experimental data (Delmas, 2007). The fitted values are the following: 0.16 for the sill λ 1 , 0.08 for the nugget effect λ 2 and 12 for the range λ 3 . Some perturbations will then be applied to theses reference values.  Figure 7.

THE DRAW DISPERSION STUDY
The variograms parameters used to fit the sample variograms have been introduced in Section 3.
We also defined the domain of variation of the variogram parameters. Our aim is now to quantify the influence of theses variations on the storage.
First the variation of the response due to the different draws with the same variogram parameters, called draw dispersion, is examined. Indeed the quantity of stored carbon is influenced by the initial configuration of the reservoir, which differs from one run to another. The induced variability can be interpreted as a repeatability error. It is then summed up by a variance term, noted σ 2 geostatistical . In practice, a hundred simulations of the same random field have been generated, representing 100 different examples of the same reservoir with, for the silicates: the range λ 1 is equal to 0.16, the nugget effect λ 2 to 0.08 and the sill λ 3 to 12. For the porosity, we use the variogram: (5) and the permeability is deduced using the linear modal of corregionalisation between the gaussian transform of porosity and permeability field (Diedro, 2009). The reactive transport software, run for each of the initial fields, give 100 different output values, of which the variance is evaluated. The results for solid stored carbon and mean porosity are presented here after (Fig. 8-11).

Solid Stored Carbon
The effect of the draw dispersion on the stored carbon (normalised to a volume of 1 m 3 of the mesh) in the solid matrix of the reservoir is presented in Figure 9. The difference between the simulations increases during the first thousand years and then stays constant. At time t = 2 000 years, the mean m is 22 196 moles and the standard deviation σ is 1 297.3 moles, which corresponds to a coefficient of variation of σ/m about 6%.

Mean Porosity
It can be seen in Figure 10 that the porosity firstly decreases; the reactions of precipitation dominate and lead to some reservoir occlusion. After 300 years the inverse trend is observed; the porosity begins to grow with dominant dissolution reactions.
The variance observed on curves can be evaluated and is illustrated in Figure 11. At time t = 2 000 years, the standard deviation σ is equal to 1.65 × 10 -5 and the mean m to 3.29 × 10 -4 . That corresponds to a coefficient of variation of 5.6%.
In this preliminary study, the quantification of the impact of the draw dispersion shows a light effect on carbon storage and porosity. In the next section, the draw dispersion is compared to the dispersion induced by the variation of the variogram parameters.

Choice of the Design
To detect the most influential parameters among the 8 selected parameters, i.e. the parameters that are the most linked to the storage, the variation of the quantity of carbon stored is quantified when the factor varies from its minimal value to its maximal value (the range of variation was defined in Sect. 3). However, each parameter x considered separately can be non-influent on the storage whereas they can influence the output if they are considered simultaneously. Therefore, the chosen model that links the inputs to the output is the following (p = 8): The coefficients x i represent the different parameters of the variogram (Tab. 1).
The coefficients correspond to half the expected variation of the answer when the factor varies from its minimal value to its maximal value (Helbert, 2006;Linder, 2005). All the factors are preliminary transformed to have the same scale of variation [-1, 1]. The value -1 corresponds to the minimal value of the initial scale and the value +1 to its maximal. The coefficients β, ..., β p are expressed in the unit of output and can be compared to each other.
The best designs of experiments to estimate that kind of models are full factorial designs that include all possible combinations of high and low levels, i.e. 2 8 experiments. These designs allow the estimation of the main effects and also the estimation of the effects of all interactions (product of two, three, …, eight factors). Moreover, the effects are estimated independently from each other and with the minimum of variance (Montgomery, 2005). However, the reactive transport software is time consuming and the realisation of the full factorial design cannot be considered. The chosen design, presented in Table 1, is a fractional factorial design of resolution IV. It contains a part of the experiments of the full design. Here, only 32 runs are launched among the 256 possible runs. These runs are chosen so that the estimation of the main effects (coefficients β, ..., β p ) are still independent. However, the consequence of reducing the number of runs is an augmentation of the estimation variance and that several effects are confounded. For example, in Table 1, the column 6, that contains the levels low and high of the sill of silicates, is confounded to the product of the first 4 columns. Thus, if we detect that this column is influent on the storage, we won't be able to distinguish if it comes from the effect of the sill or from the effect of the product x 1 x 2 x 3 x 4 . However, in most of the time, the effect of a product of more than 3 factors is neglectable. That is why the present design allows the 1-11 1-1-11 1 1-11-11111 1-11-1-11-1-1 1-1 -11111-1 1-1-11-11-11 1-1-1-11-1-11 1 -1-1-1-1-1 1 -1 -11111-11-1 -1 1 1 1 -1-1-1 1 -11-11-11 1 1 -11-1-11-11 1 -1 1 -1-1-1-1-1-1 estimation of the main effects (effects of principal factors) that are confounded with products of more than 3 factors but not of all the effects of interactions (themselves confounded with other interactions). The help of physical sense is then crucial to detect what effect takes the highest part of the estimated confounded effect. More precisely, the following relations can be obtained between interactions: x 1 * x 3 = x 7 * x 8 (7) x 1 * x 5 = x 6 * x 7 (8) x 1 * x 7 = x 3 * x 8 = x 5 * x 6 (9) x 1 * x 8 = x 3 * x 7 (10) The influence of interactions can't be mathematically separated. However, some of the interactions are naturally neglected exploiting knowledge based on the physical phenomena. In this study, all the interactions situated on the right hand of the above equations are neglected (Montgomery, 2005).
Some effect or several interactions are confounded as we have shown in Equations (7-11) below. The effect of the sill, the range and the nugget effect of the petrophysic (porosity and permeability) and the sill and the spot of the dolomite field could clearly be determined.

Analysis of the Variance
One simulation has been run for each of the 32 experiments. The two operational responses (porosity and stored carbon) are observed. In this section, our aim is to analyse the variance of the 32 observations in order to detect which of the 8 variographic parameters influence most the storage.
However, the seed of the simulation has not been fixed during the study. Thus, the total observed variance is not only due to the effect of variographic parameters variations but also to the draw dispersion. Before interpreting the results, we compare the total variance including both effects (variables and draw dispersion) to the variance induced by the draw dispersion only.
This approach makes the strong assumption that the draw dispersion does not depend on the variographic parameters. The variance above in Section 3 is the same for each run of the design. However, this is not the case since the considered distributions are lognormal and the draw has a stronger effect when the sill is high.

Results of the Analysis
As for stored carbon dioxide is concerned, the total variance is equal to 4.44 × 10 9 whereas the variance due to the draw is 1.68 × 10 6 . This indicates that more than 99% of the variance is due to variation of variogram parameters and only 1% to the repeatability error (draw dispersion).
The situation is similar for the average porosity where: σ 2 totale = 7.11 × 10 -8 (12) σ 2 geostatistical = 2.72 × 10 -10 (13) σ 2 variables = 7.09 × 10 -8 These results suggest that the noise induced by draw dispersion is very low so that the impact of variogram parameters can be estimated properly. Thanks to the results obtained on the 32 simulations, the following models for the total quantity of stored solid carbon and for the average porosity are fitted: where A is the set of the chosen interactions, where the coefficient β i quantifies half of the mean output variation when x i goes from its minimal to maximal value and where ε is a white noise.

Estimation of the Model for Stored Carbon
The evolution with time of the various estimated coefficients is observed in Figure 12. The value of each coefficient determines the impact of the corresponding parameter on the stored carbon dioxide. The major role of three of them, x 5 and the interaction x 4 x 5 , is clearly noticed. They correspond to the parameters controlling the size of the spots of dolomite within the rock. Therefore, the threshold and the sill of the dolomite are the parameters the most significant on the storage of carbon dioxide in solid form in the reservoir.
The other parameters have a negligible effect. At t equal 2 000 years, the following relation between stored carbon and variogram parameters is the following: NC =59766+37616 x 5 -35 434 x 4 -35 559 x 4 x 5 (16) The number of moles of stored carbon dioxide (NC) increases with x 5 , i.e. the size of the spot of dolomite but decreases with x 4 , the sill of dolomite.

Physical Interpretation
The stored carbon dioxide in the matrix is distributed in the form of carbonate precipitation, following the equation: determined by (Lasaga, 1998). The precipitation pattern is proportional to its reactive surface S m (Brosse et al., 2005). The greater the initial reactive surface is, the higher the dolomite precipitation will be. This confirms the importance of the spots of the dolomite within the rock for the geological storage of carbon dioxide in the matrix. An increasing size of the dolomite spots (with lower mean concentration within) indicates that the mineral offers an improved contact area with all other minerals. This promotes the reactions in the reservoir.

Results on the Average Porosity
Similarly to the previous response, the spot and the sill of the dolomite are the parameters that most influence the variation of the average porosity in the reservoir. The coefficients have a similar but inverse pattern compared to the carbon dioxide storage (see Fig. 13). The estimated equation is: Time evolution of the various coefficients of the input parameters and their interaction: X 1 , X 2 , X 3 , …, X 4 : X 8 having an impact on the carbon dioxide stored.

Figure 13
Time evolution of the various coefficients of the initial parameters and their interaction: X 1 , X 2 , X 3 , …, X 4 : X 8 having an impact on the variation of the average porosity in the reservoir.
The variation of the porosity increases with the sill of the dolomite field X 4 and decreases with the size of the spots. The stored carbon dioxide and the precipitated mineral reduce inside porosity. An increase of the dolomite spot areas results in a decrease of the average porosity in the reservoir. We also notice the opposite effect for sill: the higher the sill is, the higher the average porosity level becomes. From this experimental design study, the importance of the dolomite spots is shown. The last section of the paper deals with additional experiments carried out with constant dolomite content.

Simulation with Constant Quantity of Dolomite
To complete and better understand the results, a series of experiments with constant quantity of dolomite is performed. We simulate a dolomite field with a mean equal to 0.025, a standard deviation equal to 0.05 and a variogram range of 4 meshes. We then apply three different thresholds, a high one at 0.0267, a medium one at 0.025 and a low one at 0.024. This produces three different sizes of the dolomite spots. The fields are then normalized so that the global quantity of dolomite is maintained constant between three runs; therefore, the surface of the spots is different (Fig. 14). We obtain three different cases with various surface spot of dolomite. This is explained by the imposed constraint of constant dolomite. The transition from small to larger spots leads to a decrease of the average proportion of dolomite in the spots from 25% to 3%. In the 3 examples below section (5.4), the neutral aquifer meets the reservoir containing free CO 2 at residual saturation in conditions such as the aquifer is the only mobile phase.  Carbon stored (in mol) for three different sizes of the spots of dolomite within the rock: small spots (CS_low_dol), medium spots (CS_mean_dol) and large spots (CS_High_dol).
An important parameter is the geometry of the dolomite spots. This result (Fig. 15) is in agreement with the overall kinetics of mineral precipitation, as proportional to the mineral surface (Lasaga, 1998); carbon storage is thus proportional to a surface parameter. It should be noted however that the surface parameter is defined at the scale of several decimetres or metres and not defined at the scale of individual mineral grains, the size and surface of which are kept constant in the computer code.

CONCLUSION
The main objective of this study was to show the influence of the initial variability of the geological reservoirs on the underground storage of the CO 2 . We completed a study of experimental design by varying initial parameters around their mean value to determine their impact on two operational parameters such as the quantity of stored carbon in the reservoir and the average porosity. An important effect of the size of the dolomite spots within the rock on the random value of carbon has been identified. Whatever the value of the sill of the variogram of the dolomite field, the quantity of stored carbon increases when the range or the extent of the dolomite spot, increases. As specified in (Lasaga, 1998), the precipitation of a mineral is linked to its specific surface. A great extent of the spot of dolomite facilitates an important reactive surface and leads to more precipitation of carbonates in the reservoir.
We can also notice a similar effect obtained on the variation of the average porosity. The precipitation of a mineral reduces the average porosity. Results presented in this paper show that the average porosity decreases when the extent of the dolomite spot (i.e. the range of the threshold field) increases. It is important to notice that, with the exception of the variations of the extent of dolomite spots, small errors in the modeling of the variogram for parameters the nugget effect, the sill and the range of petrophysics and the range and nugget effect of silicates (i.e. X 1 to X 3 and X 6 to X 8 ) do not cause significant variations on the stored solid carbon and change of the average porosity. Two additional studies have been performed to complete results obtained by before. The first study concluded that the variability inferred by the statistical dispersion is weaker than the variability inferred by perturbations on the variogram parameters. The second study consisted in experimenting with a constant quantity of dolomite and confirms the important role of the spot of dolomite. All the results show that a great extent of the spot of dolomite leads to an important reactive surface, which facilitates the precipitation reactions in the reservoir.
Additional petrophysical and mineralogical data would be necessary to better calibrate our geostatistical models. This would allow the direct computation of horizontal variograms and thus avoid the use of a multiplicative coefficient for the transition from the vertical information to that for the horizontal fields. Therefore, we are confronted to the lack of data due to the difficulties to collect the mineralogical data. The results of the experimental design highlighted the important role of the spots of dolomite within the rock on the storage of carbon. So, we could try to better understand the important role of the dolomite spots using physical considerations by assessing the reactive surface associated to a given size of the spots as a function of the grain size. For an equal quantity of mineral, a slow moving fluid can quickly react when compared to a fast moving fluid. These remarks could encourage us for the future studies to define a non dimensional number using the range of the variability, the chemical kinetic and the speed of the fluid. For further studies, it will also be benefit to compare the incertitude due to the choice of chemical system.
The model proposed in the paper has been restricted to the influence of the nature and proportions of minerals on some CO 2 storage parameters; it remained conceptual in so far as it is too simple to be readily applied. However, it may be the basis for future studies where the uncertainty on the proportions of minerals will be combined with other types of uncertainties so as to be more realistic.