Thomas Kremer
email: kremer.tho@gmail.com.

Cristian Vieira

Alexis Maineult

ERT monitoring of gas injection into water saturated sands: modelling and inversion of cross-hole laboratory data

Keywords: CO 2 and N 2 injection, electrical resistivity tomography, saturated sands, transport modeling, cellular automaton, electrode configuration

HAL is a

INTRODUCTION

Since the beginning of the industrial era, anthropogenic activities such as combustion of fossil resources and deforestation have resulted in a continuous increase of atmospheric CO 2 concentrations [START_REF] Neftel | Evidence from polar ice cores for the increase in atmospheric CO2 in the past two centuries[END_REF][START_REF] Keeling | Geoelectrical methods for monitoring geological CO2 storage: first results from cross-hole and surface-downhole measurements from the CO2 SINK test site at Ketzin (Germany)[END_REF], which have recently passed over the symbolic threshold of 400 ppm [START_REF] Showstack | Carbon dioxide tops 400 ppm at Mauna Loa, Hawaii[END_REF]. It is very likely that these elevated concentrations participate to major environmental impacts such as climate change [START_REF] Manabe | Sensitivity of a global climate model to an increase of CO2 concentration in the atmosphere[END_REF][START_REF] Cox | Acceleration of global warming due to carbon-cycle feedbacks in a coupled climate model[END_REF] and ocean acidification [START_REF] Doney | Ocean acidification: the other CO2 problem[END_REF]. Carbon capture and storage (CCS), which consists of capturing the CO 2 emitted by high emissions industrial units and storing it in deep geological formations, is considered today as a viable strategy to mitigate these impacts.

To ensure the long term efficiency of the process and prevent environmental impacts associated with CO 2 leakage into nearby aquifers, it is vital to develop and utilize monitoring tools that allow the detection and quantification of any CO 2 transfer out of the targeted geological formation. Several approaches can be employed to achieve this goal. Some of them focus on the monitoring of the CO 2 plume itself (generally in the supercritical state), and aim to detect any change occurring in this area. Other strategies are directed toward the monitoring of the geological formations located above the plume up to the surface, in which case the goal is to detect a prospective CO 2 transfer within these structures (which would likely happen under the gaseous state if the pressure and temperature conditions allow it). The present study falls in the second category. Several approaches can be employed; in particular, geophysical methods can be used to study in an indirect way various physical properties of the subsurface, and detect possible variations of these properties due to the presence of leaking CO 2 . Among them, the electrical resistivity tomography (ERT) method is of particular interest because the electrical resistivity of geologic materials is sensitive to the nature, the composition and the distribution of fluids in the porous space. Furthermore, the costs of implementation are significantly lower than those of seismic methods (generally used by the oil and gas industry), thus providing a monitoring strategy with a higher temporal discretization.

To date, the ERT method has been successfully applied by numerous authors for CO 2 detection into water saturated geological media [START_REF] Xue | Estimation of CO2 saturation from time-lapse CO2 well logging in an onshore aquifer, Nagaoka, Japan[END_REF][START_REF] Bergmann | Fluid injection monitoring using electrical resistivity tomography-five years of CO2 injection at Ketzin, Germany[END_REF][START_REF] Auken | Imaging subsurface migration of dissolved CO 2 in a shallow aquifer using 3-D time-lapse electrical resistivity tomography[END_REF]. Depending on the physico-chemical state of the studied medium, two main processes can affect the ERT response of the ground. In some cases, the gaseous CO 2 is observed to cause a global increase of the medium resistivity, as it acts as an electrically insulating phase within the porous space. Such effects have been observed and well documented through field experiments [START_REF] Xue | Estimation of CO2 saturation from time-lapse CO2 well logging in an onshore aquifer, Nagaoka, Japan[END_REF]Kiessling et al., 2010;[START_REF] Bergmann | Fluid injection monitoring using electrical resistivity tomography-five years of CO2 injection at Ketzin, Germany[END_REF][START_REF] Carrigan | Electrical resistance tomographic monitoring of CO2 movement in deep geologic reservoirs[END_REF] and at the laboratory scale [START_REF] Nakatsuka | Experimental study on CO2 monitoring and quantification of stored CO2 in saline formations using resistivity measurements[END_REF][START_REF] Breen | Bench-scale experiments to evaluate electrical resistivity tomography as a monitoring tool for geologic CO2 sequestration[END_REF][START_REF] Bosch | Core-scale electrical resistivity tomography (ERT) monitoring of CO 2-brine mixture in Fontainebleau sandstone[END_REF]. In other cases, it has been observed that the bulk resistivity of the ground can decrease when important dissolution processes are occurring [START_REF] Peter | Investigation of the geochemical impact of CO2 on shallow groundwater: design and implementation of a CO2 injection test in Northeast Germany[END_REF][START_REF] Vialle | Percolation of CO2 rich fluids in a limestone sample : Evolution of hydraulic, electrical, chemical, and structural properties[END_REF][START_REF] Auken | Imaging subsurface migration of dissolved CO 2 in a shallow aquifer using 3-D time-lapse electrical resistivity tomography[END_REF]. The general sequence of events starts with the dissolution of some of the gaseous CO 2 into the saturating water, which becomes more acidic, thus enhancing the potential dissolution of minerals. All the newly dissolved ionic species lead to a more electrically conductive water and hence to a less resistive medium. In some cases, dissolution and desaturation effects can be observed simultaneously [START_REF] Xue | Estimation of CO2 saturation from time-lapse CO2 well logging in an onshore aquifer, Nagaoka, Japan[END_REF][START_REF] Barrio | CO 2 Migration Monitoring Methodology in the Shallow Subsurface: Lessons Learned from the CO 2 FIELDLAB Project[END_REF][START_REF] Yang | Monitoring CO 2 migration in a shallow sand aquifer using 3D crosshole electrical resistivity tomography[END_REF].

In this paper, we present metric scale experiments which aim to simulate a gaseous CO 2 leakage through a relatively simple medium (well characterized sand saturated with water) and monitor the resistivity variations induced by the gas circulation. Although strong differences obviously exist, these experiments can be seen as an attempt to approach the conditions of a borehole-to-borehole ERT monitoring network through which a CO 2 leakage path would form, in geological conditions close to those of a saturated aquifer with a relatively homogeneous structure. Two different sands were used in order to study in parallel carbonated environments (in which strong dissolution processes can occur) and silica based environments (in which dissolution processes are very weak). The experiments were performed using gaseous CO 2 and gaseous N 2 in order to study separately the "insulating phase" effect, when gas invades the porous space and no dissolution processes occur. Before the experiments, a preliminary study was conducted to investigate different ERT measuring protocols and choose an optimal one. Then, we present the results of the laboratory experiments. In view of the results, a model was designed under the form of a cellular automaton with stochastic behavior, in order to simulate gas injection and circulation through the medium. Coupling this model with the software R3t [START_REF] Binley | R3t manual[END_REF], we were able to reproduce the ERT responses observed during the experiments and to infer the geometry of the gas plume, which changes depending on the grain size distribution of the sand. Finally, in a third phase, the cellular automaton model was used as a tool to study numerically the sensitivity of the electrodes configurations and to develop new, more efficient acquisition protocols.

MATERIALS AND METHODS

Experimental set-up

To perform the experiments, an 80-cm-high, 48-cm-diameter cylindrical Plexiglas tank was used (Fig. 1a). Electrodes were arranged in 10 rings of 6 cylindrical electrodes fixed on plastic rods attached on the inner wall of the tank. The vertical ring spacing was 5 cm, and the lowest ring was located 18 cm above the bottom of the tank (Fig. 1b). The electrodes were not placed exactly on the wall of the tank due to the diameter of the rods (1 cm) and the space between each rod and the internal wall (3 mm). The electrodes (4 mm in diameter and 5 mm in length) were made of stainless steel to prevent from corrosion. Gas injection was performed using a valve placed at the bottom center of the tank. The gas was injected through a cylindrical air diffuser (10 cm in height and 12 cm in diameter, Fig. 1c) placed on the axis of the tank, between the planes 10 and 8. A small pipe connects the diffuser to a valve placed at the bottom center of the tank, and linked to gas reserves. All these components were made of PVC to minimize electrical disturbances.

Porous medium

Experiments were performed with two types of sand. The first one is a silica sand from Fontainebleau (France) and the second is a carbonate sand obtained from commercial aggregates. The porosity values of the sands were estimated using decimetric columns filled with sand, by measuring the mass difference between the unsaturated state and the water saturated state. The permeability estimates were obtained by applying different hydraulic heads values on the saturated columns and monitoring the hydraulic flow. The structural and hydraulic characteristics of the sands are listed in Table 1, where errors were calculated as the standard deviation of the results obtained on 4 independent measurements. To create a homogeneous distribution in the tank and avoid bubble trapping, the sand was slowly, regularly poured into water using a sieve whose mesh size is slightly larger than the largest grain size. Sand and water were then added and stirred simultaneously so as to always keep a free water layer of a few centimeters thick at the top of the medium, which was removed once the filling procedure was completed. Such procedure ensured that no significant air bubbles remained compared to the volume of the tank and the amount of gas injected during the experiments. The air diffuser was placed after the lowest ring of electrodes was reached. We avoided the formation of unfilled cavities in the zones between the vertical plastic rods and the internal wall of the tank, since they could act as preferential gas migration pathways. No compaction process was applied during the filling procedure except in these zones.

Gas injection and dissolution processes

We employed dinitrogen (N 2 ) and carbon dioxide (CO 2 ) during the experiments.

Dinitrogen (N 2 ) presents a very low solubility into water (about 0.016 volume per volume at a temperature T = 0° C), and does not react with the solid phase, while carbon dioxide (CO 2 ) has a solubility into water at 0° C of about 0.88 and can react with the solid matrix once dissolved in water. In the conditions of the experiments, these gases can be considered as perfect insulators (their resistivity is about 15 orders of magnitude higher than the resistivity of the water used during the experiments). It is important to note that we worked with an opened system, where the injected gas is allowed to pass through and flow out of the medium (the tank is not covered up). Gas injection is done according to the following protocol: a first injection phase is performed at a rate of about 200 L•h -1 for at least 30 min, creating a steady state with multiple gas pockets spread in the tank and gas migration pathways that seem stable, as indicated by the observation of bubbles appearing continuously at the same locations at the surface of the medium. Then, the injection rate is increased to 400 L•h -1 . Once a steady-state is reached (about 30 min later) the ERT acquisition starts and the injection is sustained until the whole acquisition sequence is performed, which lasted about 1.5 hours.

Experiments involving N 2 injection allowed studying the case where changes of resistivity are solely induced by an increase of the gas saturation, and where no dissolution processes occur. On the other hand, during CO 2 injection, chemical reactions due to carbonic acid might occur depending upon mineralogy [START_REF] Xu | Mineral sequestration of carbon dioxide in a sandstone-shale system[END_REF][START_REF] Busch | Effects of physical sorption and chemical reactions of CO 2 in shaly caprocks[END_REF][START_REF] Raistrick | Carbon dioxide-water-silicate mineral reactions enhance CO 2 storage; evidence from produced fluid measurements and geochemical modeling at the IEA Weyburn-Midale Project[END_REF]. When put in contact with water, part of the gaseous CO 2 dissolves and reacts with water forming carbonic acid, which will itself dissociates, providing additional ionic species into the water. Those processes can be summarized by the following equations.

In the case of a solid phase composed of calcium carbonate, the dissolved CO 2 will also enhance the dissolution of part of the solid matrix, further increasing the ionic concentrations in the saturating water.

CO 2 and minerals dissolution can hence generate an input of dissolved species which modify the water resistivity value. Therefore, in experiments involving CO 2 injection, two phenomena can possibly affect the effective resistivity of the media: 1) water saturation level variation due to gas invasion, which would provoke an increase of the effective media resistivity and 2) dissolution processes, which would result in a decrease of the media resistivity. In the following, we will refer to these two phenomena respectively as the "insulating phase" effect and the "dissolution processes" effect.

Data acquisition and processing

The acquisition device was a Terrameter SAS 1000 connected with a 64-electrode multiplexer ES-1064 by ABEM. Injected current consists of square waves of alternating polarity of 20 mA to avoid polarization and associated noise [START_REF] Parasnis | Principles of Applied Geophysics[END_REF]. The acquisition delay was set to 0.1 s and the acquisition time to 0.2 s.

To perform a tomography reconstruction, we measured the potential differences resulting from current injection for a large set of electrodes arrays. Finally, to obtain an estimation of the resistivity distributions from the measured resistances, inversions were carried out with the software R3t [START_REF] Binley | R3t manual[END_REF], as it is suitable to the medium geometry and experimental conditions. R3t is based on a regularized objective function with weighted least square. Neumann conditions are specified at the boundaries where no current flow is allowed (see [START_REF] Binley | R3t manual[END_REF] for more details). The solution of the inverse problem is obtained through a Gauss-Newton procedure, as described in [START_REF] Binley | DC resistivity and induced polarization methods[END_REF]. In this study we used the background regularization scheme described in [START_REF] Binley | R3t manual[END_REF], where the regularization parameter is optimized each iteration by carrying out a line search. The convergence of the inversion process occurred only when the root-mean-square (RMS) error, as defined in [START_REF] Binley | R3t manual[END_REF], was below 1.

The electrodes were modeled as point sources. This was considered an acceptable hypothesis given the results from [START_REF] Rücker | The simulation of finite ERT electrodes using the complete electrode model[END_REF]. According to their study, the errors related to the point source assumption become important when the ratio between the electrode spacing and electrode length overcomes 0.2. In the present study this ratio never overcomes 0.1 for all measurements.

As explained in [START_REF] Binley | Regularised Image Reconstruction of Noisy Electrical resistance Tomography Data[END_REF]Binley et al., ( , 1996a[START_REF] Binley | Flow pathways in porous media: electrical resistance tomography and dye staining image verification[END_REF], [START_REF] Binley | DC resistivity and induced polarization methods[END_REF] and, [START_REF] Daily | Electrical resistance tomography[END_REF], the characterization of the error is required for ERT modeling to assess its effect during reconstruction of the true resistivity distribution. The error can be inferred from repeatability and reciprocity. Repeatability is related to the precision at which measurements are obtained under identic conditions. Reciprocity is based on the principle that switching positions between injection electrodes A and B and measuring electrodes M and N should yield the same resistance value.

The measurement error was estimated using both repeatability and reciprocity principles [START_REF] Binley | DC resistivity and induced polarization methods[END_REF][START_REF] Daily | Electrical resistance tomography[END_REF] during a preliminary experiment (N 2 injection into saturated siliceous sand). One estimation was made before starting the gas injection and another one while gas injection was occurring but when gas migration pathways seemed stable. The error inferred from reciprocity was found to be always less than 1% in both cases. This value was considered low enough to be ignored and reciprocity was therefore not assessed during the experiments presented in this paper. Only the repeatability error was taken in account in the inversion process.

Optimal acquisition arrays

An important issue in ERT surveys relates to the sensitivity of the electrode configurations used for the acquisition and, by extension, to the sensitivity of the measuring protocol. The optimization of the electrode configurations on the surface of the soil (surface ERT) has been studied by several authors [START_REF] Dahlin | A numerical comparison of 2D resistivity imaging with 10 electrode arrays[END_REF][START_REF] Christensen | Monitoring CO2 injection with cross-hole electrical resistivity tomography[END_REF][START_REF] Szalai | Depth of investigation and vertical resolution of surface geoelectric arrays[END_REF]. To our knowledge, most cross-well ERT monitoring strategies implement dipole-dipole or bipole-bipole arrays for acquisition [START_REF] Pezard | Time-lapse downhole electrical resistivity monitoring of subsurface CO 2 storage at the Maguelone shallow experimental site (Languedoc, France)[END_REF][START_REF] Commer | Time-lapse 3-D electrical resistance tomography inversion for crosswell monitoring of dissolved and supercritical CO 2 flow at two field sites: Escatawpa and Cranfield, Mississippi, USA[END_REF][START_REF] Yang | Monitoring CO 2 migration in a shallow sand aquifer using 3D crosshole electrical resistivity tomography[END_REF], Preignitz et al., 2013) sometimes mixing both arrays [START_REF] Hermans | Quantitative temperature monitoring of a heat tracing experiment using cross-borehole ERT[END_REF]. A review of pros and cons of some classical cross-borehole arrays are also presented in [START_REF] Bing | Cross-hole resistivity tomography using different electrode configurations[END_REF]. In the case of measurements at the surface of a cylinder, [START_REF] Binley | DC resistivity and induced polarization methods[END_REF] recommend to use all linearly independent combinations. We chose to investigate this issue for our particular case, by testing several protocols on simple media with As a first test, for each protocol, the resistances measured on a medium composed of silica sand fully saturated with water were compared with the resistances values obtained through forward modelling of the same medium. The silica sand was saturated with a NaCl solution concentrated at 0.56 g.L -1 . The resistivity of such medium was measured thanks to a test cell, which yielded a value of 38 Ω•m. This resistivity value was then used to perform a forward modeling of the acquisition. Fig. 3 displays the results obtained with the different protocols tested. Concordance between measured and modeled resistances is much better when the cross-rods protocol is used than when the other "classical" configurations are used (bipole-bipole, dipole-dipole or Schlumberger). Given these results, we used exclusively the cross-rods protocol in the experiments described hereafter.

EXPERIMENTAL RESULTS

As explained in the previous section, two types of gas (N 2 and CO 2 ) and two types of sand (silica and carbonate) were used in this study, which resulted in the realization of four different experiments, one for each gas/sand combination. In the following, we will refer to each experiment using the initial of the sand type and the gas type formula. For example, C/N 2 will be referring to the experiment involving the injection of dinitrogen into saturated carbonate sand.

ERT measurements on pure water

Gas injection into carbonate sand

Fig. 5 shows the resistivity distribution obtained for the experiments involving carbonate sand (C/N 2 and C/CO 2 ). For these experiments, the carbonate sand has been saturated with a NaCl solution whose initial concentration was 0.25 g.L -1 ( 18.5 Ω.m).

During and after sand saturation, the resistivity of the water was found to decrease slowly, until reaching a stable value of about 14 Ω‧ m at 17.5°C, which is most likely due to the dissolution of some carbonated particles during and immediately after the saturation process.

Fig. 5a shows the resistivity distribution of the initial state, before gas injection, which is quite homogeneous except for some small artifacts in the vicinity of the electrodes. A slightly higher resistivity zone can also be observed in the central part of the tank, which might be a consequence of the plastic gas diffuser presence.

Fig. 5b and 5c show the results of experiments C/N 2 and C/CO 2 , respectively. For both experiments, the resistivity distributions obtained are no longer uniform. The N 2 injection results in the formation of a high resistivity anomaly, likely to be induced by gas accumulation into the porous space. This anomaly exhibits a u-shape stretching from the gas diffuser to the surface of the tank and represents a significant resistivity increase of about 50 % compared to the initial state. The injection of gaseous CO 2 results, on the contrary, in a decrease of the resistivity values compared to the initial state. This anomaly is also u-shaped, centered with the axis of the tank, and stretching from the top of the diffuser to the surface of the sand pack. It is likely that the injection of gaseous CO 2 into a reactive medium such as carbonate sand leads to strong dissolution processes which increase the electrical conductivity of the saturating water and hence reduce the medium resistivity. In this case, even if gas accumulation is also occurring, the impact of dissolution processes is stronger and dominates the resistivity response of the medium.

We recall that the gas injection sequence consists of a first 30 min duration injection at a rate of 200 L‧ h -1 followed by a second 30 min duration injection at a rate of 400 L‧ h -1 . At that point the "steady-state" has been reached (gas circulation patterns seem stable through visual inspection of the surface of the medium) and the ERT acquisition starts and lasts for about 1.5 hours. We note that the "steady-state" assumption is valid for experiment C/N 2 but not necessarily for experiment C/CO 2 , where dissolution processes might still be occurring while the ERT acquisition is still in progress, thus leading to a kind of average resistivity distribution over the duration of the ERT acquisition (even though most dissolution processes are likely to occur before the beginning of ERT acquisition, given the fast kinetics of these reactions). We applied another inversion method in order to better extract the changes due to gas injection only, and remove the effects induced by the presence of the gas diffuser. This method is called the ratio-inversion. As explained in [START_REF] Miller | Application of time-lapse ERT imaging to watershed characterization[END_REF] and [START_REF] Daily | Electrical resistance tomography[END_REF], a common approach for detecting changes from a background image is to normalize the data according to the following equation:

where is the post-injection data, is the zero state data, and is the data obtained through forward modeling of a homogeneous medium. The inversion process is then applied using this normalized data set (in our case a set of resistances values). Fig. 6 show the results obtained with the ratio inversion method for experiment C/N 2 and C/CO 2 . The shapes of the anomalies are globally very similar to those obtained with a "simple" inversion, i.e u-shaped, centered with respect to the tank's vertical axis and stretching from the top of the diffuser to the surface of the sand pack. Based on these results, the ratio inversion method does not seem to provide significantly more accurate results, since the same shapes can be observed with a simple inversion. However, if one would intent to use these results in order to quantify the quantity of N 2 or CO 2 injected, it should be done with the ratio inversion results, which are supposed to take in account only the resistivity structures induced by gas injection, and exclude those present in the background resistivity distribution (like the gas diffuser).

Gas injection into silica sand

Fig. 7 shows the resistivity distribution obtained for the experiments involving silica sand (S/N 2 and S/CO 2 ). For these experiments, the sand was saturated with an electrolyte whose resistivity was 21 Ω‧ m at 18.5°C, which resulted in a higher background resistivity than for the carbonate experiments. As for the carbonate sand experiments, the initial state resistivity distribution (not shown here) was again quite homogeneous, exhibiting values mainly between 40 and 50 Ω‧ m and only few artifacts in the vicinity of the electrodes. Fig 7a

shows the resistivity results obtained with a simple inversion. The injection of gaseous N 2 induces again the formation of a high resistivity anomaly, with values above 70 Ω‧ m, that we can attribute to the invasion of the porous space by the insulating gaseous phase. However, the shape of the anomaly is quite different since it does not seem to reach the upper part of the tank, as if most of the gas was being trapped in the lower part. The ratio-inversion method (Fig. 7b) manages to reduce some part of the anomaly that was likely due to the presence of the diffuser, indicating that gas accumulation in the lowest part of the tank might not be as significant as the simple inversion results might suggest. However, the shape of the anomaly remain mainly horizontal, not stretching toward the top of the sand pack and hence being quite different from the results obtained with carbonate sand. Fig. 7c shows the results of the ratio inversion for experiment S/CO 2 . Again, we observe a high resistivity anomaly, with values even higher than for experiment S/N 2 . Even though it is likely that some of the CO 2 has been dissolving into the water, the "insulating phase" effect seems to be dominating in this case. In this case, the zone around the diffuser still exhibits high resistivity values after ratio inversion, which might be attributed to significant gas accumulation around the diffuser itself. Also, the anomaly again seems restricted to the lower part of the tank, showing horizontal components rather than the u-shape observed with carbonate sand. These different shapes suggest that the gas circulation patterns might be different depending on the type of sand used for the experiment.

Summary

These results confirm that the ERT method is a suitable monitoring technique to detect gas injection into geological media. Yet, the shape and the magnitude of the resistivity anomalies seem to depend on the nature of the sand, the type of gas injected and probably on the initial electrical conductivity of the saturating water. When the silica sand is involved, the anomaly seems to be restricted to the lower part of the tank, as if gas could not reach or accumulate in the upper volume. When the carbonate sand is used, the anomaly exhibits a ushape centered with the axis of the tank, and spreading from the top of the diffuser to the surface of the tank. Also, for a highly reactive system such as carbonate sand and gaseous In the following discussion, we want to address two main questions raised in view of these experimental results. 1) In the case of slightly reactive systems, why does the increase in resistivity exhibits different anomaly shapes depending on the type of sand involved. 2) Is it possible to develop new measuring protocols that will further improve the ability of the method to capture the shape and magnitude of these anomalies?

Note that these two questions seem to overlook the case where dissolution processes dominate the electrical response of the system whereas it is relevant in many scenarios of CO 2 leakage in the subsurface. We justify this choice because these particular cases have important temporal and kinetic aspects that are not addressed nor captured by the experiments presented previously, which rather focus on steady-state situations. However, they can both still bring pertinent answers to some of the possible CO 2 leakage scenarios, as well as to any other problems involving gas circulation in the subsurface or which demand the detection and delimitation of a resistivity anomaly.

In the next section, we will address question 1) and discuss how the sand grain size and the gas injection rate might be important factors controlling the geometry of the gas circulation plume, which results in varying resistivity anomaly shapes. We present a numerical model that was developed to challenge this hypothesis. It is structured as a stochastic cellular automaton and allows to simulating the injection of gas into saturated sands and computing the resulting electrical resistivity distribution of the medium. Finally, in section 5, we use this model to address question 2) and investigate the development of new measuring protocols that would improve the detection of these anomalies.

MODELLING

4.1 The effect of grain size and gas injection rate on gas plume geometry.

Apart from issues related to the geological storage of CO 2 , the way that gas spreads into a saturated geological medium has also been studied during the development of the "air sparging" method [START_REF] Böhler | Air injection and soil air extraction as a combined method for cleaning contaminated sites-observations from test sites in sediments and solid rocks[END_REF][START_REF] Johnson | An overview of in situ air sparging[END_REF], a remediation technique which consists of injecting air into or below a shallow contaminated area. While the air flows towards the surface, the contact with the pollutant triggers biodegradation processes and volatilization, allowing some of the pollutant to flow toward the surface under gaseous form, along with the injected air. The context of air-sparging is relevant here because it involves the circulation of a nearly non-reactive gas (air), which is why many studies dealt with laboratory experiments very similar to those described in this study, especially those involving N 2 injection.

At the end of the 1990s, most of the authors working on the air-sparging method agreed that, in a saturated medium, the air circulates according to two main patterns. The "bubble flow" consists of a discontinuous circulation under the form of bubbles [START_REF] Wehrle | In-situ cleaning of CHC contaminated sites : model-scale experiments using the air injection (in-situ stripping) method in granular soils[END_REF][START_REF] Sellers | Air sparging model for predicting groundwater cleanup rate[END_REF]; while the "channel flow" describes a circulation through continuous pathways, often called channels, which evolve independently from each other and are separated by zones fully saturated with water [START_REF] Clayton | A field and laboratory investigation of air fingering during air sparging[END_REF][START_REF] Elder | Air channel formation, size, spacing, and tortuosity during air sparging[END_REF].

The laboratory studies from [START_REF] Ji | Laboratory study of air sparging : Air flow visualization[END_REF], [START_REF] Baker | Review of factors affecting in situ air sparging[END_REF] and [START_REF] Brooks | Fundamental changes in situ air sparging patterns[END_REF] have shown that the main discriminating factor between these two flow modes is the grain size of the porous medium. [START_REF] Brooks | Fundamental changes in situ air sparging patterns[END_REF] gathered from several studies the different kind of flow that can be observed depending on the mean grain size of the studied medium.

This comparison (see their figure 1) suggests that for a grain size below 1 mm, the air circulation follows the channel flow type, and for a grain size above 2 mm only the bubble flow is observed. In between, a transition behavior is observed which consists of discontinuous channels, and that later authors have been calling "mixed flow" or "slug flow".

For thin sands (d < 0.21 mm), [START_REF] Peterson | Air-flow geometry in air sparging of fine-grained sands[END_REF] identify a third type of gas flow that they call "chamber flow". It involves the existence of continuous channels, which will accumulate in some location of the porous medium, following mainly horizontal components, hence forming gas accumulation "chambers" connected between each other through continuous channels that are essentially vertically oriented. [START_REF] Brooks | Fundamental changes in situ air sparging patterns[END_REF] have also studied the influence of the gas injection rate upon the type of flow encountered. Their observations suggest that low injection rates favor the "slug flow" or "mixed flow" type, whereas high injection rates allow the gas to form continuous channels. The different kind of flow observed depending on grain size and injection rate are summarized in their table 1. Regarding the experiments described in this paper, the silica sand and the carbonate sand maximum grain size are respectively 400 µm and 2 mm, and the injection rate is high with respect to the references from [START_REF] Brooks | Fundamental changes in situ air sparging patterns[END_REF]. We can therefore safely assume that the gas circulation during the experiments followed a "channel flow" type or at the very least a "slug flow", and this assumption has been considered throughout the following modeling work.

Modelling gas injection into saturated sands

The model developed in this paper has been built as a cellular automaton using Matlab. A 3D parallelepiped matrix is generated, whose dimension is , being the diameter of the cylindrical tank used during the experiments (Fig. 8). With this model, each term of the 3D matrix corresponds to a 1 cm 3 volume of the tank that we refer to as a cell. The model consists of simulating gas circulation patterns through the cells, then defining an effective water saturation value for each cell and finally estimating the associated resistivity value using Archie's law [START_REF] Archie | The electrical resistivity log as an aid in determining some reservoir characteristics[END_REF].

Modelling gas circulation patterns

As explained previously, based the grain size of the sand used during the experiments, it can be considered that gas circulation follows the "channel flow" type. The first step of the modeling process is hence to define the number of gas channels created during the injection, their dimensions and their respective trajectory. [START_REF] Selker | The local geometry of gas injection into saturated homogeneous porous media[END_REF] describe the formation of gas channels as a competition between restriction forces imposed by the capillary pressure and buoyant forces produced by the density difference between gas and water. Assuming that channels are cylindrical and that their extremity is hemispheric, [START_REF] Selker | The local geometry of gas injection into saturated homogeneous porous media[END_REF] define the channels radius value to be where is the interfacial tension between gas and water (N/m), Δρ is the density difference between gas and water, and g is the gravity acceleration (m.s -2 ). The computation of R in the case of N 2 into water saturated sands yields a value of approximately 0.5 cm. This implies that the channel diameter scale is similar to the scale of one cell of our cellular automaton grid. [START_REF] Selker | The local geometry of gas injection into saturated homogeneous porous media[END_REF] also propose an estimate of the number of gas channels created, that depends on the injection rate and on the physical and hydraulic properties of the gas and granular medium, respectively:

where Q is the gas injection rate (m 3 .s -1 ), μ is the dynamic viscosity of the gas (kg.m -1 .s -1 ) and κ is the permeability to gas (µm 2 ) of the medium, that we estimate using the permeability to water values given in Table 1. Equation 9 yields about 300 channels for the silica sand and about 100 channels for the carbonate sand.

Stochastic evolution of channel trajectory

Once the number and dimensions of gas channels have been estimated, the next step is to define the geometry of the gas plume, i.e the trajectory of the gas channels. Looking back at Fig. 5b, the resistivity anomaly suggests that the gas spreads following a parabolic shape or "u-shape", which has indeed been observed in several studies [START_REF] Ji | Laboratory study of air sparging : Air flow visualization[END_REF][START_REF] Semer | TECHNICAL NOTE an experimental investigation of air flow patterns in saturated soils during air sparging[END_REF][START_REF] Selker | The local geometry of gas injection into saturated homogeneous porous media[END_REF]. Fig. 7, however, suggests that the gas is mainly concentrated in the lower part of the tank. This kind of behavior has been less frequently reported, but we assume in the following that it can be related to the "chamber flow" described by [START_REF] Peterson | Air-flow geometry in air sparging of fine-grained sands[END_REF]. Indeed, the chamber flow pattern has been reported for fine grain sands (d < 210 µm),

which is approximately consistent with the grain size of the silica sand used in the experiment (210 µm < d < 400 µm). It is therefore possible that the gas accumulation in the lower part of the tank consists of one "chamber" spreading laterally, and that gas is flowing toward the surface through vertical channels whose concentration is not high enough to produce significant resistivity anomalies. This hypothesis is supported by visual observations that were made during similar gas injection experiments using the same silica sand (see Kremer et al., 2016a and their Figure. 9b).

Our first objective in this modeling work is to implement the parabolic shape gas circulation mode and the "chamber" flow mode, to see if they can respectively explain the resistivity anomalies observed during experiments C/N 2 and S/N 2 . For [START_REF] Selker | The local geometry of gas injection into saturated homogeneous porous media[END_REF], the parabolic circulation pattern can be explained by the fact that in a relatively isotropic medium, the probability that a gas channel moves toward one lateral direction or another is the same. Therefore, the probability distribution associated with the position of a gas channel in the horizontal plane follows a Gaussian law, and the channel concentration is maximal at the injection point.

To model this behavior, a stochastic approach has been adopted. It consists of simulating the trajectory of each gas channel independently as the successive invasion of neighboring cells. Each channel begins at a cell chosen randomly just above the gas diffuser.

The trajectory of the channel is then constrained to 2 possibilities: 1) move up in one of the cells located above which share at least one edge with the present cell (9 possibilities); 2) move laterally in one of the cells sharing at least on edge with the present cell (8 possibilities).

The first probabilistic event "gas moves vertically" is named V, of probability P V , while the second probabilistic event "gas moves horizontally" is named H, of probability P H . As the universe is restricted to these two events we have For a given value of P H (or P V ), a selection with equal probability is performed to determine whether the channel moves vertically or horizontally. Then another selection is done to determine which of the 8 (or 9) cells will be invaded. This is repeated until a cell at the top of the tank is reached. After that, we consider that the channel is stable and that gas is flowing out of the tank. This stochastic process is applied independently to each gas channel created.

At the end of the simulation, the number of channels passing through each cell is recorded.

The shape of the gas plume is then controlled by the value of P H , as shown on Fig. 9 where, as expected, high P H values result in a wider gas plume, with a higher number of channels in the middle part.

This way of simulating gas channel evolution preserves the Gaussian distribution of pathways described by [START_REF] Selker | The local geometry of gas injection into saturated homogeneous porous media[END_REF], as shown in Fig. 10a. To model the chamber flow accumulation describe by [START_REF] Peterson | Air-flow geometry in air sparging of fine-grained sands[END_REF], the same stochastic approach is used, except that the probability P H is no longer fixed, but is a function of the vertical position of the current cell.

where is the actual position of the channel extremity, is the diffuser vertical location and denotes the height of the surface. Such a definition implies that the probability of horizontal movement is high close to the diffuser location and decreases exponentially as the channel evolves toward the surface. It results in channel accumulation in the lower part of the tank, while the channel density is lower in the upper part (see Fig. 10b and10c), thus reproducing the concept of "chamber flow accumulation".

Defining saturation level and resistivity values

As explained previously, given the definition of gas channel dimensions (equation 6), the channel diameter is approximately of the same order as the model cell size. Therefore, we consider in the following that several channels passing through the same cell only result in one channel passing through the cell. A fixed gas saturation value Sg is then attributed to all cells invaded by gas. This value is chosen randomly between 0.4 and 0.5, as observed by [START_REF] Chen | Computed tomography imaging of air sparging in porous media[END_REF]. The saturation values distribution is then smoothed by averaging the value of each cell with the value of its closest neighbors. Finally, the resistivity value of each cell is calculated using Archie's law [START_REF] Archie | The electrical resistivity log as an aid in determining some reservoir characteristics[END_REF]. The porosity values of the sands are known, and the cementation exponent was chosen to be 1.3 (unconsolidated sands). Note that surface conduction is not taken in account with Archie's law, which is acceptable given the conductivity value of the fluid which is in the range where surface conduction is negligible compared to bulk conduction [START_REF] Glover | A generalized Archie's law for n phases[END_REF]. Fig. 11 displays the vertical sections of the saturation values obtained with the model and the associated resistivity distribution for the "parabolic" circulation case and for the "chamber accumulation" case. The gas circulation patterns are modeled assuming that they follow a "channel flow" type, as discussed previously. We parametrize their geometry and density based on the study conducted by [START_REF] Selker | The local geometry of gas injection into saturated homogeneous porous media[END_REF] which describes the formation of gas plumes into saturated sands. The way these channels circulate within the tank is controlled by probabilistic laws, described in appendix A. They are key parameters that allow modelling different gas plume geometries and reproduce those observed during the experiments. In order to model the injection of N 2 into the carbonate sand, we use probabilistic parameters that produce a u-shape (or parabolic) gas plume (Fig. 11a), as suggested by the experimental results (Fig. 5 and6), and as observed in several studies about gas injection into sands with similar characteristics [START_REF] Ji | Laboratory study of air sparging : Air flow visualization[END_REF][START_REF] Semer | TECHNICAL NOTE an experimental investigation of air flow patterns in saturated soils during air sparging[END_REF][START_REF] Selker | The local geometry of gas injection into saturated homogeneous porous media[END_REF]. In the case of N 2 injection into silica sands, the resistivity shape observed in the experiments reminds the "chamber flow" geometry observed by [START_REF] Peterson | Air-flow geometry in air sparging of fine-grained sands[END_REF] on similar fine grain sands. The probabilistic parameters are therefore adapted to reproduce this specific geometry (Fig. 11d). Then, gas saturation values are assigned based on [START_REF] Chen | Computed tomography imaging of air sparging in porous media[END_REF], smoothed by averaging with the neighboring cells (Fig, 11b and11e), and resistivity values are finally computed using Archie's law (more detail in appendix A). The results are 3D representations of the resistivity distribution within the tank, as shown on Fig. 11c and11f.

Comparing model and experimental results.

The software R3t was then used to: 1) calculate the forward model associated with these two synthetic resistivity distributions, considering the same electrodes locations and measuring protocol (cross-rods) as during the experiments; 2) calculate the inverse model using the so computed synthetic resistances. Fig. 12 compares the inversion results obtained based on the synthetic models (b and d) with those obtained based on experimental measurements (a and c). The resistivity range is very similar in both cases (between 60 and 90 Ω•m), and the shapes of the high resistivity anomalies also match, thus confirming that in silica sands, the gas is likely to be circulating following the chamber accumulation flow type whereas in carbonate sands, the parabolic gas plume is a plausible description of the gas circulation patterns. These results suggest that, in saturated sands, the grain size distribution is likely to be the main factor controlling the shape of the gas circulation patterns, and hence of the resistivity anomaly that is observed. The differences between the observed shapes are quite significant. Grain size distributions that result in high permeability values (e.g. 16.5 µm 2 for carbonate sand) will favor the occurrence of u-shaped gas plumes that have a strong vertical component. For lower permeability values (e.g. 4,5 µm 2 for silica sand), the gas plume takes the form of an accumulation chamber with dominating horizontal components. One could therefore expect that grain size variation within an aquifer will produce specific resistivity anomaly shapes alternatively dominated by vertical and horizontal orientations. These observations highlight the fact that knowing the grain size distribution of sandy aquifers is very important if one wants to anticipate the shape of the gas circulation patterns, and hence of the prospective resistivity anomaly that would be observed through ERT monitoring of the aquifer.

DISCUSSION

In light of these results, we investigated the possibility to develop new measuring protocols that would improve the capacity of the ERT method to detect and capture the shape and magnitude of the resistivity anomalies produced by gas circulation into saturated sands.

Indeed, even if some information is available to help anticipate the shape of the expected resistivity anomaly, it is possible that different cross-hole protocols would not be equally efficient into detecting and retrieving the gas plume geometry because of their different sensitivity characteristics.

As a reminder, the cross-rods protocol was initially chosen to perform the experiments because it was proven to be more efficient than the bipole-bipole, dipole-dipole and Schlumberger protocols (see section 2). Here we attempt to design even more efficient protocols than the cross-rods. To do so, our approach consists of analyzing the sensitivity distribution of all the electrode configurations that composed the cross-rods protocol, so as to combine them in a new, more optimal organization. The respective efficiency of these new protocols is then tested using the cellular automaton model, by simulating rough synthetic resistivity distributions as well as realistic gas circulation patterns.

Electrode configurations families.

The cross-rods protocol consists of performing all possible electrode combinations for current injection and voltage measurement, in each of the 10 planes (rings) of electrodes instrumenting the tank. Electrodes from different planes are never used together for a single measurement. This protocol hence results in several electrode geometric combinations. To be more specific, as a plane contains 6 electrodes, the number of combinations possible for 4 electrodes measurements yields 45 possibilities (not including reciprocals). Among these 45 possibilities, one can distinguish 8 different geometric families, which are represented in Fig. 13 and that we name "classes" 1 to 8.

Sensitivity distribution of the different families.

Given the different geometries of these electrodes combinations, they present different sensitivity distributions that we chose to estimate using the following procedure. The objective is to perform forward modeling using the cross-rods protocol on two cylindrical media (representing the cylindrical tank) to which a different resistivity distribution is applied.

We attributed to the first medium a homogeneous resistivity distribution of 25 Ω•m (same order of magnitude than the saturated sand). The second medium has the same background These sensitivity distributions were estimated for each electrode geometric configurations and regrouped in Fig. 14, where high relative variations values denote a higher sensitivity area.

The resolution of this mapping is relatively poor (resolution equivalent to the cylindrical intrusion diameter), but was very quick to implement, and nevertheless yields quantifiable and valuable information regarding the sensitivity behavior of the different electrodes configurations used in the cross-rods protocol.

Fig. 14 shows that the geometries can be separated into two different groups. Geometries n°1, 4, 6 and 7 have a relatively high sensitivity distribution, rather homogenously distributed and symmetrical with respect to the center of the tank (excepting n°4). Among them, geometry n°6 is of particular interest since it presents the highest sensitivity values at the center of the tank, and high values (> 10%) on a diametrically symmetrical fringe. Another group composed of geometries n°2, 3, 5 and 8, rather show high sensitivity values at the lateral areas of the tank. In particular, geometry n°3 has the highest values distributed over approximately a third of the tank circular surface. This suggests that, taken together, the 3 electrode arrays that possess geometry n°3 provide a high sensitivity scan of the entire tank surface, except in the center.

Another electrode configuration which deserves special attention is geometry n°8. This specific set-up yields extremely high relative variations compared to the others (in the order of 10 ^4%). This is due to the numerical nature of the study, and is directly related to the fact that the injection and measuring dipoles are perpendicular. Because of this, in a perfectly homogeneous medium, perpendicular dipoles imply that the current lines and the equipotential lines are perpendicular, yielding a null potential difference and therefore a null resistance. In our case, the obtained resistance is not exactly 0 because of the irregularities of the mesh (the tank is not perfectly cylindrical), but remains very small (of the order of 10 -3 Ohm). Then, including a different resistivity object in the environment still leads to the measurement of small resistances, in the order of 1 Ohm, which results in a relative variation of several order of magnitude. Even though these considerations depend on the numerical nature of this study, it has also been observed in the real case that the relative variations induced by gas injection in a relatively homogeneous medium are higher for the measurements associated with this perpendicular electrode geometry. In the following sections, we investigated further this specific configuration (and others), to verify if this specific behavior could play an important role in the reconstruction of resistivity contrasts and anomalies.

Designing new protocols

In this section we present several numerical tests performed in order to explore the behavior of some specific electrodes geometries and assess their influence on the quality of the inversion results. We also aim at creating new protocols that may possess an improved sensitivity and resolution compared to the cross-rods protocol. To do so, we performed a series of tests, which consisted of applying different protocols for the computation of the forward and then inverse problem, on specific resistivity distributions composed of a high resistivity object (100 Ω•m) included into a lower resistivity background (25 Ω•m) and comparing the quality of the inversion problem results between each protocol. The high resistivity object was chosen as a cross, composed of two perpendicular parallelepiped volumes (Fig. 15), so as to assess the capacity of the tested protocol to retrieve both horizontal and vertical structures.

Test 1: influence of each electrode geometry on the cross-rods protocol performance In order to investigate the importance of each of the 8 electrode geometries on the global performance of the cross-rods protocol, we performed a first test consisting of removing all measurements associated to a specific electrode geometry and compare the ability of these "truncated" protocols to retrieve the cross shaped object. The results are presented on Fig. 16.

As expected, removing some measurements from the protocol deteriorates the quality of the inversion results. However, the degree of deterioration varies depending on which electrode geometry is excluded. Geometries n° 1, 2 and 3 seem to play an important role since their removal results in a strong weakening of the results quality compared to the full cross-rods protocol. The worst result seems to be obtained when geometry n°3 is absent from the protocol. Geometry n°6 may be the most important one, since without the associated measurement the inversion process cannot converge.

Test 2: measurements with electrodes in different planes.

The purpose of this test is to analyze the contribution of measurements involving electrodes from different planes. Given the results from the previous test, we chose to develop these inter-planes protocols using electrode configuration geometries n°3 and n°6. We define an inter-plane measurement as an electrode configuration where one current injection electrode and one potential measurement electrode are located in a different plane from the two other electrodes. Fig. 17a illustrates this definition. It shows that one electrode geometry can be derived in 9 additional measurements (without considering reciprocals). A protocol was designed which included all the inter-planes measurements possible using geometries n°3 and n°6. However this protocol gathered 900 measurements, when the initial cross-rods protocol has only 450 measurements. In order to have a comparable number of measurements and thus truly compare the global efficiency of the protocol, we designed a split inter-planes protocol, for which each initial geometry is not derived in 9 additional measurements but in 4, depending on the plane position of the initial measurement (Fig. 17b). By doing so the split inter-planes protocol is composed of 576 measurements, which is closer to the number of protocols composing the cross-rods protocol.

Fig. 18 compares the results of the "cross shape" test for the cross-rods, the interplanes, and the split inter-planes protocols. We clearly observe that the inter-plane protocol gives the best results. The resistivity values are higher where the cross shape is supposed to be reconstructed, and the horizontal part of the shape is better retrieved with inter-planes protocols than with the cross-rods one. For a smaller amount of measurements (split interplane protocol), the results are slightly worse but still remain much better than the cross-rods protocol.

Test 3: random construction of protocols

Based on these results, we worked on building protocols using measurements from the crossrods protocol and from the inter-planes protocols associated with geometry 6 and 8. As it is quite difficult to identify which specific configuration brings the most information, we chose as a first step to build protocols with 450 measurements by randomly selecting 150 electrodes configurations in the cross-rods and in the 6 and 8 inter-planes protocols. We generated 10 different protocols and compared their results to the cross-shape test. Only 3 over 10 produced better results than the cross-rods protocol, among which one showed particularly good results that we named Random. By studying which electrode geometries were picked-up initially, we noticed that these 3 protocols have in common a high amount of measurements from geometry 2, 3, 4, 5, 6 and 8, and a low amount from geometry 1 and 7. Considering that the random pick-up from 6 and 8 inter-planes protocol already bring a sufficient amount of geometries 6 and 8, we chose to restrict our random pick up in the cross-rods protocol to geometries 2, 3, 4 and 5, in proportions of 30, 60, 30 and 30 measurements, respectively.

Among the 10 protocols generated with this method, 8 demonstrated a superior behavior at the cross-shape test than the cross-rods protocol, and one called the Random Constrained was selected as the most efficient.

Test on new geometries and realistic resistivity distribution

From the three tests previously described, three protocols emerged which showed a better capacity than the cross-rods protocol to reconstruct the cross-shape object: the split interplanes, the Random, and the Random constrained protocols. Fig. 19 compares the ability of these protocols to reconstruct other geometries such as a non-centered cross (Fig. 19a) and a diagonal rod (Fig. 19b). In these two cases, the three new protocols still produce a better restitution of the initial resistivity distribution, and the highest resistivity values are reached with protocol Random constrained. Fig. 19c and 19d compare the results of forward modeling and inversion processes performed on resistivity distributions that represent the "chamber flow" and "u-shape" gas plumes simulated with the cellular automaton model. In these cases, results are not anymore so consistent. The Random constrained protocol still demonstrates a satisfying behavior (2 nd best for both distributions), but is not anymore the most efficient. For the chamber flow distribution, the Random protocol is superior, whereas for the u-shape distribution the cross-rods protocol is the most efficient protocol. These results demonstrate that the efficiency of the protocol to detect a gaseous plume greatly depend on the geometry of the plume. Here we found that Random constrained seems to be quite polyvalent whereas

Random and cross-rods are respectively the best choices to image a "chamber flow" or a "ushape" geometry.

CONCLUSION

The laboratory experiments confirmed the ability of the ERT method to detect, locate and image CO 2 and N 2 plumes into water saturated sands. The experiments validated the occurrence of two different effects leading to opposite resistivity variations compared to the background distribution: the "insulating phase" effect and the "dissolution processes" effect.

Also, the resistivity anomalies produced by N 2 injection exhibit different shapes depending on the type of sand used during the experiments. Parabolic shapes are observed with carbonate sands, whereas more flattened shapes are observed with silica sand.

Using a gas circulation model designed as a cellular automaton with stochastic behavior, we managed to reproduce the resistivity distributions observed during the experiments involving N 2 injection and hence infer the geometry of the gas circulation plumes. The interpretation of these results with respect to existing literature on two-phase flow into granular media strongly support the hypothesis that the particle size distribution of the sand is an important factor controlling the geometry of gas circulation plumes, and hence of the resistivity anomalies obtained during the ERT monitoring. Sands with larger grain diameter ( ) will tend to produce u-shape or parabolic anomalies, whereas in smaller grain distributions ( ), gas can accumulate within horizontal chambers which will produce flattened anomalies.

Finally, we developed new acquisition protocols and produced a variety of tests to study their characteristics and performances. These investigations demonstrated that the capacity of one acquisition protocol to retrieve the true resistivity distribution after data inversion can vary significantly depending on the shape of the resistivity anomaly that is to be detected. Therefore, in order to optimize the measuring arrays for the monitoring of CO 2 or any other gas transfer through the subsurface, one must carefully assess the prospective geometry of the gas circulation plume, and hence must gather relevant information on the geological structures composing the subsurface.

Fig. 1 .

 1 Fig. 1. Experimental set-up. Picture of the cylindrical tank equipped with the electrodes (a) and of the gas diffuser (c). Sketch of the tank and position of the electrodes (b).

  uniform resistivity distribution. The arrays bipole-bipole, dipole-dipole, Schlumberger and the protocol called cross-rods configuration (recommended by[START_REF] Binley | DC resistivity and induced polarization methods[END_REF] were tested. For the bipole-bipole array (Fig.2a), current adjacent electrodes are located in the same plane whilst potential electrodes are located on the same position in the upper or lower plane. In dipole-dipole (Fig.2b) the current electrodes are placed in the same rod and the potential electrodes in the adjacent rod, but the separation between A and B and M and N covered all possible combinations. In the Schlumberger array (Fig.2c) all the electrodes were placed in the same rod where the distance between A and M equals the distance between N and B. The cross-rods configuration consists of all possible 4-electrodes combinations located in the same plane (on the same ring), for each plane. All the possible 4-electrodes combinations are represented on Fig.13. The main advantage of this protocol is that it forces the current lines to pass through the sand pack, enhancing the penetration depth.

Fig. 2 .

 2 Fig. 2. Scheme of the "classical" acquisition arrays: bipole-bipole (a), dipole-dipole (b) and Schlumberger (c). Letters A and B denote the positions of the current injection electrodes and M and N that of the potential measurement electrodes.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Comparison between measured and modeled resistances. (a,b): Schlumberger protocol; (c,d): bipole-bipole; (e,f): dipole-dipole; (g,h): cross-rods configuration. Left column: measured (orange) and modeled (black) resistances values for each measurement. Right column: modeled vs measured resistances.

Fig. 4

 4 Fig.4shows the resistivity distribution obtained after inversion when applying the

Fig. 5 .

 5 Fig. 5. Resistivity distribution obtained after inversion of ERT data acquired on carbonate sand saturated with water of resistivity 14 Ω.m: at the initial state (left), during N 2 injection of (middle), and during CO 2 injection (right). The black rectangle indicates the position of the diffuser.

Fig. 6 .

 6 Fig. 6. Resistivity distribution obtained after application of ratio inversion on the data from experiments C/N 2 (left) and C/CO 2 (right). The resistivity anomalies remain u-shaped, stretching from the top of the diffuser to the top of the tank, and are centred with respect to the axis of the tank.

Fig. 7 .

 7 Fig. 7. Inverted resistivity distribution for gas injection in saturated silica sand. a: N 2 injection, simple inversion. b: N 2 injection, ratio inversion. c: CO 2 injection, ratio inversion. The position of the air diffuser is given by the open rectangle.

  CO 2 (C/CO 2 ) the dominant physico-chemical processes are linked to the dissolution of the gas and of the solid matrix, leading to a decrease of the resistivity, as observed for example by[START_REF] Peter | Investigation of the geochemical impact of CO2 on shallow groundwater: design and implementation of a CO2 injection test in Northeast Germany[END_REF] or[START_REF] Vialle | Percolation of CO2 rich fluids in a limestone sample : Evolution of hydraulic, electrical, chemical, and structural properties[END_REF]. Yet, for non-reactive or slightly reactive systems (S/N 2 , S/CO 2, C/N 2 ), dissolution processes are weak and the dominant mechanism is the desaturation of the porous space, which increases the resistivity of the medium , as observed for example in the studies conducted by[START_REF] Nakatsuka | Experimental study on CO2 monitoring and quantification of stored CO2 in saline formations using resistivity measurements[END_REF],[START_REF] Breen | Bench-scale experiments to evaluate electrical resistivity tomography as a monitoring tool for geologic CO2 sequestration[END_REF] or[START_REF] Bosch | Core-scale electrical resistivity tomography (ERT) monitoring of CO 2-brine mixture in Fontainebleau sandstone[END_REF].

Fig. 8 .

 8 Fig. 8. Isometric view of the parallelepiped matrix. Green cells have been assigned the value 1 and represent the studied medium, while transparent cells have been assigned the value NaN and correspond to the outside of the tank.

Fig. 9 .

 9 Fig. 9. Vertical cross-sections of the parabolic shape circulation model, for several values of probability P H. The color code represents the total number of channels passing through each cell.

Fig. 10

 10 Fig. 10. a) Histogram representing the lateral distribution of channels in the parabolic shape model. The distribution is globally Gaussian. b) Vertical cross-section of the accumulation chamber model, and c) corresponding vertical distribution of channels.

Fig. 11 .

 11 Fig. 11. Vertical and horizontal cross-sections representing the distribution of channels, saturation values and resistivity values, respectively, for the parabolic circulation model (a,b,c) and the chamber accumulation model (d,e,f).

Fig. 12 .

 12 Fig. 12. Comparison between the inversion results obtained from the experimental data and from the cellular model, for the accumulation chamber (a, b) and the parabolic cases (c, d), respectively.

Fig. 13 .

 13 Fig. 13. Schematic representation of the different four electrodes arrays families encountered in the cross-rods protocol.

Fig. 14 .

 14 Fig. 14. Sensitivity distributions obtained for each class of electrode array geometry. Letters C-, C+ and P+-, P+, denote the current injection electrodes and potential measurement electrodes locations, respectively.

Fig. 15 .

 15 Fig. 15. Representation of the cross-shape resistivity intrusion embedded into the homogeneous resistivity cylindrical medium.

Fig. 18 .

 18 Fig. 18. Results from test 2. Comparison between the cross-shape test using a) the cross-rods protocol (450 measurements), b) the inter-planes protocol (900 measurements) and c) the split inter-planes protocol (576 measurements).

Table 1 . Characteristics of the sands

 1

	Sand
	Silica	Carbonate

ACKOWLEDGEMENTS

We gratefully acknowledge the support of the CO 2 research program co-funded by ADEME, TOTAL and SCHLUMBERGER and their permission to publish this paper, and duely note that the views expressed here are those of the authors. We also thank Dr. Preignitz, Dr G. Donno, as well as an anonymous reviewer, whose comments greatly help to improve this article.