Mohammad Ali Iravani

Hossein Davarzani
email: h.davarzani@brgm.fr

Jacques Deparis

Stéfan Colombano

Nicolas Philippe

Benjarese Oniangue

Roger Guérin

Alexis Maineult

Hossein Davarzani 1'

Experimental study of electrical complex resistivity in a 2D multiphase porous medium under non-isothermal conditions: application to soil remediation monitoring

Keywords: Induced polarization, modeling, pollution, resistivity, environmental

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Pollution of soil due to dense non-aqueous phase liquid (DNAPL) spilling has become a crucial issue of concern during the last decade, causing numerous geo-environmental (e.g., mitigation of DNAPLs on aquifers) problems [START_REF] Habiyakare | The implementation of genetic algorithm for the identification of DNAPL sources[END_REF]. In polluted sites, DNAPLs sink below the water table and cause dissolved plumes in the subsurface [START_REF] Pankow | Dense chlorinated solvents and other DNAPLs in groundwater: History, behavior, and remediation[END_REF][START_REF] Power | Evaluating four-dimensional time-lapse electrical resistivity tomography for monitoring DNAPL source zone remediation[END_REF]. After many years of research, soil contamination due to DNAPLs is still one of the main challenges for the environmental and geophysical communities [START_REF] Kavanaugh | The DNAPL remediation challenge: Is there a case for source depletion? Environmental protection agency Washington DC[END_REF]. DNAPLs, including coal tar (CT), contaminate soils and water widely when released in the subsurface, mainly from industrial parks. DNAPLs are part of a group of immiscible fluids that pass through the saturated and unsaturated zones due to gravity and capillary forces, via migration passages made of pores, until stopping at an impermeable stratum [START_REF] Illangasekare | Experimental study of movement and distribution of dense organic contaminants in heterogeneous aquifers[END_REF][START_REF] Power | Evaluating four-dimensional time-lapse electrical resistivity tomography for monitoring DNAPL source zone remediation[END_REF].

Various remediation techniques can be used in these sites to depollute them. These include pump and treat, thermal, bioremediation, and chemical methods [START_REF] Beyke | In situ thermal remediation of DNAPL and LNAPL using electrical resistance heating[END_REF][START_REF] Hadley | Groundwater remediation: The next 30 years[END_REF]. Pump-and-treat is the common and first remediation technique to remove the free product. After pumping, the residual saturation can be removed using biochemical remediation techniques.

Geophysics has shown the value of the long-term monitoring of clean-up processes (e.g., [START_REF] Brewster | Observed migration of a controlled DNAPL release by geophysical methods[END_REF][START_REF] Wilson | DNAPL mapping by ground penetrating radar examined via numerical simulation[END_REF]. Of all the geophysical methods, electromagnetic methods and especially electrical methods are useful for studying the processes involved during clean-up and degradation of pollutants [START_REF] Chambers | Noninvasive monitoring of DNAPL migration through a saturated porous medium using electrical impedance tomography[END_REF][START_REF] Mao | Geophysical monitoring of hydrocarbon-contaminated soils remediated with a bioelectrochemical system[END_REF][START_REF] Deng | Application of spectral induced polarization for characterizing surfactant-enhanced DNAPL remediation in laboratory column experiments[END_REF], 3D imaging of a jet fuel spill [START_REF] Orozco | Delineation of hydrocarbon contaminants with multi-frequency complex conductivity imaging[END_REF][START_REF] Helene | Identification of leachate infiltration and its flow pathway in landfill by means of electrical resistivity tomography (ERT)[END_REF]. The electro-geophysical methods, especially direct current (DC) resistivity, induced polarization (IP) techniques, and spectral induced polarization (SIP) have proved their potential in geological conditions, in saturated conditions (e.g., [START_REF] Weller | Estimating permeability of sandstone samples by nuclear magnetic resonance and spectral-induced polarization[END_REF][START_REF] Slater | Electrical-hydraulic relationships observed for unconsolidated sediments in the presence of a cobble framework[END_REF][START_REF] Revil | Predicting permeability from the characteristic relaxation time and intrinsic formation factor of complex conductivity spectra[END_REF]Kremer et al., 2016b;[START_REF] Sparrenbom | Investigation of chlorinated solvent pollution with resistivity and induced polarization[END_REF], unsaturated conditions (e.g., [START_REF] Maineult | Variations of petrophysical properties and spectral induced polarization in response to drainage and imbibition: a study on a correlated random tube network[END_REF], oil, gas, and mining industries (e.g., [START_REF] Friedrichs | Controlled Source Electro-Magnetic Mapping development and evaluation of its application to mineral exploration[END_REF], soil pollution and remediation monitoring (e.g., [START_REF] Pearce | Electrical properties of soils at ism radio frequencies, from 1 to 40 MHz[END_REF]Ling and Zhang, 2017;Liao et al., 2018) and civil engineering (e.g., [START_REF] Susanto | Mechanical, electrical and microstructural properties of cement-based materials in conditions of stray current flow[END_REF]. In a multiphase porous medium, various combinations of different parameters such as temperature, saturation, level of compaction, and pore fluid properties can have a similar effect on the measurements (e.g., [START_REF] Sen | Influence of temperature on electrical conductivity on shaly sands[END_REF][START_REF] Denicol | Effects of water salinity, saturation and clay content on the complex resistivity of sandstone samples[END_REF]. Therefore, studying the effect of each parameter to determine their role and ideally, their related simultaneous effects on geophysical parameters is necessary. [START_REF] Sadowski | Clay-organic interactions[END_REF] performed the first laboratory study on different clay samples with different pollutants and different concentrations. He found that the SIP method could be used to identify organic contaminants in clayey soils for different concentrations. [START_REF] Börner | Contamination indications derived from electrical properties in the low frequency range 1[END_REF] contributed to laboratory research on the effects of organic liquid contaminations on different rock types and sandstones. They demonstrated that both the real and imaginary parts of complex resistivity are affected by the properties of the contaminant filling the pores. They noted that in some cases the various contaminants could be differentiated using the imaginary part and emphasized the effect of interactions between pore fluid characteristics and pores' internal surface structure.

The existence of organic contamination can affect capacitive conduction in the IP response due to the surface reaction between contamination and soil grains [START_REF] Vaudelet | Mapping of contaminant plumes with geoelectrical methods. A case study in urban context[END_REF]. IP responses in these contaminated soils are related to the saturation [START_REF] Breede | Spectral induced polarization measurements on variably saturated sand-clay mixtures[END_REF], distribution, concentration, and types of contaminants present in soil [START_REF] Breede | Spectral induced polarization measurements on variably saturated sand-clay mixtures[END_REF]. Previous studies demonstrate that the presence of contaminants (Gasoline and leachate results of acidification by organic and carbonic acids) in the soil increases resistivity [START_REF] Benson | Mapping groundwater contamination using dc resistivity and VLF geophysical methods-A case study[END_REF][START_REF] Sauck | A model for the resistivity structure of LNAPL plumes and their environs in sandy sediments[END_REF]. [START_REF] Martinho | An experimental study of organic pollutant effects on time domain induced polarization measurements[END_REF] reported an increase in chargeability due to the presence of pure organic contaminations (gasoline, isopentane, toluene, and benzene) in clayey soil. Many studies based on the SIP method have shown an increase in amplitude and phase shift due to the existence of hydrocarbons [START_REF] Vanhala | Detecting organic chemical contaminants by spectral-induced polarization method in glacial till environment[END_REF][START_REF] Börner | Contamination indications derived from electrical properties in the low frequency range 1[END_REF][START_REF] Schmutz | Influence of oil saturation upon spectral induced polarization of oil-bearing sands[END_REF]. However other SIP studies reported a decrease in amplitude and phase shift in different media consisting of different pure and dissolved hydrocarbons [START_REF] Weller | Measurements of spectral induced polarization for environmental purposes[END_REF][START_REF] Vanhala | Mapping oil-contaminated sand and till with the spectral induced polarization (SIP) method[END_REF][START_REF] Revil | Influence of oil wettability upon spectral induced polarization of oil-bearing sands[END_REF].

Recently, Iravani et al. (2020aIravani et al. ( , 2020b) evaluated how temperature and saturation changes in two DNAPLs (coal tar and chlorinated solvent) affect complex electrical resistivity and relative permittivity. They used an oven to study temperature change and a 1D drainageimbibition system to study how saturation changed these two parameters. They carried out experiments in cylindrical cells filled with glass beads (to simulate soil/porous media) and saturated with water and CT. The relationship between electrical resistivity and relative permittivity with saturation was investigated for two different coupled immiscible fluids, water/coal-tar and canola oil/salty ethanol. They found that temperature changes these two geophysical parameters linearly for the coal tar (CT) and chlorinated solvent coupled with water. The effect of saturation change was not linear, and they did state that the generalized Archie's law could fit their experimental data well. [START_REF] Colombano | Thermal and chemical enhanced recovery of heavy chlorinated organic compounds in saturated porous media: 1D cell drainage-imbibition experiments[END_REF] investigated the relationship between resistivity and saturation for chlorinated solvent. They used a drainageimbibition protocol to monitor electrical resistivity and permittivity. Their objective was to estimate residual saturation indirectly. They used the complex refractive index model (CRIM) and generalized Archie's law to model their permittivity and resistivity data, respectively.

They reported that estimated electrical resistivity presented less correlation with the experimental data compared to permittivity. They concluded that temperature and saturation affect the electrical resistivity and relative permittivity of the medium. [START_REF] Colombano | Comparison of thermal and chemical enhanced recovery of dnapl in saturated porous media: 2D tank pumping experiments and twophase flow modelling[END_REF] investigated the imbibition of a mixture of chlorinated solvents by thermal and chemical enhancements [START_REF] Colombano | Comparison of thermal and chemical enhanced recovery of dnapl in saturated porous media: 2D tank pumping experiments and twophase flow modelling[END_REF]. They found that the temperature has no effect on DNAPL recovery that could be linked to the actual properties of the oils. However, coal tar is more viscous than most chlorinated solvents. Moreover, the dynamic viscosity of coal tar decreases significantly with temperature [START_REF] Philippe | Experimental study of the temperature effect on two-phase flow properties in highly permeable porous media: Application to the remediation of dense non-aqueous phase liquids (DNAPLs) in polluted soil[END_REF]. To the best of our knowledge, no other work is reported on electrical resistivity measurements for coal tar. [START_REF] Blondel | Temporal evolution of the geoelectrical response on a hydrocarbon contaminated site[END_REF] compared lab and field results in clayey and clayey-silt soils. They quantified a resistivity increase, normalized chargeability, and decreased resistivity in measurements with oil not degraded. When it is partially degraded they saw no changes in the phase-lag and chargeability parameters. [START_REF] Gourry | Coupling of geophysical measurements and gas analyzes for the detection of pollutants organic on a former coking[END_REF] also showed that under certain conditions, oil pollution like this can increase chargeability because of high resistivity and it is rather related to wettability. This statement was not confirmed by [START_REF] Blondel | Temporal evolution of the geoelectrical response on a hydrocarbon contaminated site[END_REF]. Orlando and Palladini (2019) carried out a series of laboratory tests to study how resistivity and permittivity vary in a 2D cell in a multiphase (two-phase) system. The contamination they used was HFE-7100 (hydrofluoroether) in saturated (with pure, dissolved, and gaseous phases) grained glass bead (GB) porous media. They used geophysical measurements and physico-chemical constraints to monitor these parameters over a long period. Different geophysical techniques were used to propose empirical models to monitor water saturation in a multiphase porous media contaminated by DNAPL. Photographs were used to confirm the geophysical data. Low frequency IP can detect the pure HFE-7100 phase.

In another study, [START_REF] Orlando | Electrical permittivity and resistivity time lapses of multiphase DNAPLs in a lab test[END_REF] studied how electrical resistivity varies in ground due to the existence of DNAPL in a laboratory study. Like Orlando and Palladini (2019), they used a GB-packed medium contaminated with HFE-7100 in a saturated condition. They found that different DNAPL phases have different impacts on the physico-chemical characteristics of the contaminated medium. In their study, free and dissolved DNAPL caused a change in electrical resistivity. They concluded that electrical resistivity coupled with the GPR technique can be used for characterizing the contaminated sites because neither technique could be sensitive enough for a long time.

It is known that the presence of organic contaminants modifies the impedance of the fluidgrain interface through clay polymerization [START_REF] Olhoeft | Low-frequency electrical properties[END_REF][START_REF] Vanhala | Mapping oil-contaminated sand and till with the spectral induced polarization (SIP) method[END_REF]. The presence of organic pollutants as free phase changes the distribution of the conductive aqueous phase in the pore space (like air can do in unsaturated media). As SIP response changes with varying water/air saturation (e.g., [START_REF] Ulrich | Induced polarization measurements on unsaturated, unconsolidated sands[END_REF][START_REF] Titov | Induced polarization of unsaturated sands determined through time domain measurements[END_REF][START_REF] Binley | Relationship between spectral induced polarization and hydraulic properties of saturated and unsaturated sandstone[END_REF][START_REF] Maineult | Variations of petrophysical properties and spectral induced polarization in response to drainage and imbibition: a study on a correlated random tube network[END_REF], free-phase hydrocarbons should also produce a measurable SIP effect. A few studies have been conducted on complex resistivity laboratory measurements either on artificially contaminated soil or rock samples or on field samples from contaminated sites.

Even though the results differ in the observed effects, they proved that the SIP response of soils and rocks is measurably influenced by hydrocarbon contamination. However, a detailed understanding of how contamination affects SIP response is still lacking. Many studies did not pay sufficient attention to controlling the nature or distribution of contaminants in the porous media, and on assessing independently the effects of dissolved, adsorbed, and free-phase, leading sometimes to what appear to be contradictory results. Here, we assumed that the SIP response of the contaminant obeys a law similar to what [START_REF] Vaudelet | Mapping of contaminant plumes with geoelectrical methods. A case study in urban context[END_REF] observed above a plume of chlorinated solvents, i.e., a Cole-Cole like behavior.

The study of geophysical parameters during remediation by the thermal enhanced method is influenced by two environmental parameters: temperature due to soil heating (Iravani et al., 2020a[START_REF] Philippe | Experimental study of the temperature effect on two-phase flow properties in highly permeable porous media: Application to the remediation of dense non-aqueous phase liquids (DNAPLs) in polluted soil[END_REF] and saturation changes (Iravani et al., 2020b;[START_REF] Colombano | Thermal and chemical enhanced recovery of heavy chlorinated organic compounds in saturated porous media: 1D cell drainage-imbibition experiments[END_REF] due to contamination pumping. In previous studies, we separately investigated the impact of temperature and saturation on electrical resistivity in contaminated saturated porous media with CT (Iravani et al., 2020a(Iravani et al., , 2020b)).

In this study, we investigated the coupling of transient heat transfer in multiphase flow with saturation variation in a 2D tank. No work to date has studied transient temperature transfer and saturation changes at the same time in 2D due to its complexity. Heating in the 2D tank mimics the soil heating in the field by increasing the temperature of pollutants to decrease their viscosity. These phenomena increase pollutant mobility and so increase the pumping rate of contamination in saturated soil during depollution. Heating the samples also changes the geo-electrical parameters. The goal of the current study is to determine the saturation field in the 2D tank from resistivity measurements. All these experimental studies have been done to implement the SIP for monitoring the real remediation process of DNAPLs in the field.

It is worth noting that one of the main challenges in the monitoring of the subsurface with geophysical methods is that there are few methods (e.g., boreholes) for verifying the recorded geophysical response of the medium to validate the remediation process [START_REF] Spies | Electrical and electromagnetic borehole measurements: A review[END_REF]Keys, 2017).

Theoretical background

Spectral-induced polarization mechanisms SIP is an extension of the IP method with alternating current (AC) in the frequency domain. SIP is used to study the ratio and phase shift between voltage and current of the subsurface or a sample. Generally, electrical polarization in the range of studied frequencies (i.e., from 0.01 to 10,000 Hz) has four main mechanisms including electrode polarization, membrane polarization, Maxwell-Wagner polarization, the polarization of the diffuse layer (electrical double layer, EDL) including the Stern layer polarization (e.g., [START_REF] Breede | Spectral induced polarization measurements on variably saturated sand-clay mixtures[END_REF]. This study does not separately investigate the different mechanisms involved to characterize the effect of CT saturation and temperature changes.

Electrical complex resistivity

The electrical resistivity (ρ) in a homogeneous sample can be obtained by multiplying resistance (R) and a geometric factor (K) related to the position of the electrodes (for current injection and voltage measurement). The complex electrical resistivity reads (e.g., [START_REF] Guéguen | Introduction to the Physics of Rocks[END_REF][START_REF] Ghorbani | Contribution au développement de la résistivité complexe et à ses applications en environnement[END_REF][START_REF] Glover | Geophysical properties of the near surface Earth: Electrical properties[END_REF]: [START_REF] Olhoeft | Low-frequency electrical properties[END_REF] where ' and '' are, respectively, the real and the imaginary parts of complex resistivity and and are the amplitude and the phase shift of the resistivity response (for low phase, below a few hundred of mrads, can be simplified as . These equations can be written in terms of the electrical conductivity, (S.m -1 ) which is the inverse of the electrical resistivity.

A very first relationship to study the effect of temperature on electrical resistivity was presented by [START_REF] Dakhnov | Geophysical well logging: The application of geophysical methods; electrical well logging[END_REF] as:

(2)

where  w is the resistivity of water at temperature T and  w0 is the resistivity of water at the initial temperature T 0 , and α is the temperature coefficient of resistivity (see Iravani et al, 2020a, for instance). We used equation ( 2) to obtain α value for our experiments using thermocouples, with a reference temperature of T 0 =20°C, where T is the temperature measured by thermocouples, and  w0 is the averaged measured resistivity at 20°C. This linear relationship was defined for water but can be used for DNAPL mixtures (Iravani et al., 2020a).

In the literature, no universal model has been presented with the consideration of the physico-chemical properties of each phase in a multiphase system in regard to the frequency changes. The lack of a developed comprehensive model on this subject steers researchers to fit the experimental data of complex resistivity with the phenomenological models. One of the most used phenomenological models is presented by [START_REF] Cole | Dispersion and absorption in dielectrics I. Alternating current characteristics[END_REF] and was derived from the Debye (1929) model, which was firstly developed for relative permittivity [START_REF] Pelton | Mineral discrimination and removal of inductive coupling with multifrequency IP[END_REF]. [START_REF] Pelton | Mineral discrimination and removal of inductive coupling with multifrequency IP[END_REF] proposed a model for SIP spectra that exchanged relative permittivity with electrical complex resistivity for practical applications that describes

(3)

where  0 is the low-frequency resistivity, c is the so-called Cole-Cole exponent, m is the chargeability, and  is the mean relaxation time (m and c can be in the range of 0 and 1). If c equals 1 this equation will be similar to the Debye model. For this study, the Cole-Cole model was chosen because it has been demonstrated to be a good model to represent empirically the polarization of single grains or pores (e.g., Leroy et al. 2008;[START_REF] Niu | Connecting complex conductivity spectra to mercury porosimetry of sedimentary rocks[END_REF].

In this study for saturation prediction, the generalized Archie's law, which is developed by [START_REF] Glover | A generalized Archie's law for n phases[END_REF] based on the classical Archie's law [START_REF] Archie | The electrical resistivity log as an aid in determining some reservoir characteristics[END_REF], was used. It is defined as below Archie's law (4) and Generalized Archie's law ( 5) where (S/m) is the conductivity of the liquid inside the pores, Ф is the porosity, i is the indication of each phase and m is the cementation exponent.

In this study, we validate mixing models and calibration curves obtained from the previous study in cells and columns (Iravani et al., 2020a(Iravani et al., , 2020b) using a 2D two-phase flow and comparing the results with imaging for validation. In a previous study, Iravani et al.

(2020a) showed both resistivity and phase are related to the temperature variation in 1D columns. They found a linear relationship between temperature and resistivity. By using the Cole-Cole model they also found that temperature affects the frequency domain. Iravani et al.

(2020b) studied the saturation impact on the magnitude and phase of electrical complex resistivity. They reported the generalized Archie's law was applicable to study the multiphase medium contaminated by CT. They found that saturation change has a direct impact on electrical resistivity.

Material and methods

A laboratory-scale experiment was carried out to study temperature and saturation effects on complex electrical resistivity during the drainage and imbibition of DNAPLs in a 2D tank. The experiment was done in a saturated porous medium with an immiscible liquid pair (CT and water). In the following, we describe the essential points, experimental setup, and procedures used for this experiment.

2D tank setup

The 2D tank has dimensions of 30 cm × 50 cm × 7 cm (height, width, depth). A schematic form of the front, back, and top of the tank is shown in Figure 1. The body of the 2D tank is made of a solvent resistant material called polyvinylidene fluoride (PVDF). The front wall of the tank with a thickness of 0.8 cm is made of a transparent glass rectangle through which pictures are taken at each stage.

A stainless-steel heating element was inserted inside the left-side of the tank to produce the non-isothermal condition. Fourteen type K thermocouples (green circles) record temperatures in the tank. The thermocouples were connected to a Campbell CR-1000 acquisition system. The temperatures were measured with sampling intervals of 30 s. For the SIP measurements (Figure 1c), black and blue circles represent 30 current and 30 potential electrodes, respectively, 5 cm apart from each other. The electrical AC is injected between two stainless steel current electrodes (diameter=65 mm) made of nickel-cobalt alloy (MP35N) and the response voltage is measured between two handmade non-polarizable potential electrodes (Cu/CuSO 4 ) with a diameter of 1 cm. Non-polarizable electrodes are used in SIP studies to measure the voltage response of an injected current in a medium (e.g., Kremer et al. 2016a;[START_REF] Abdulsamad | Assessing the high frequency behavior of non-polarizable electrodes for spectral induced polarization measurements[END_REF]. In non-polarizable electrodes, current can freely pass without polarization. Among common electrode solutions, Cu/CuSO 4 was chosen based on the method proposed by [START_REF] Maineult | Electrical response of flow, diffusion, and advection in a laboratory sand box[END_REF]. The solution was made using a mixture of 72.75 g of milli-Q water (ultrapure water), 26 g of CuSO 4 , and 1.25 g of gelatin in 100 g of solution.

A heating shaker was used to mix the solution for one hour at 80°C.

The non-polarizable electrodes were placed outside the samples to avoid inflow interference that in turn can interfere with electrical current signals [START_REF] Maineult | Electrical response of flow, diffusion, and advection in a laboratory sand box[END_REF]. A porous ceramic (Ameteksi scientific instruments) was used in the interface of non-polarizable electrodes and the saturated porous media. Each non-polarizable electrode was connected to the SIP-lab IV (manufactured by Radic Research) with a shielded cable to avoid contamination of the data due to electromagnetic coupling [START_REF] Maineult | Electrical response of flow, diffusion, and advection in a laboratory sand box[END_REF]Bagotsky, 2005).

Glass beads (GB) made of silica with a diameter of 1 mm, forming a pack with a porosity of 40% ± 2%, were used as a porous medium. Different GB sizes were tested and this size was chosen to simulate a porous medium in the lab with close physical properties to field soil (e.g., the grains size and the porosity of the medium). The tap water was degassed using an ultrasonic tank (VWR Ultrasonic Cleaner -USC500D). For obtaining temperature variation effects in a non-isothermal condition, an immersion heater heats from the left side of the tank, and heat transfers to the right side. In non-isothermal conditions, there is a possibility of studying the temperature and saturation changes at the same time based on how temperature varies across the different halves of the tank and how saturation varies during imbibition from the bottom middle of the tank.

Water is injected using five pipes located at the bottom of the tank with the possibility of cutting fluid currents by using clips and clamps. The pipes are numbered and connected to a peristaltic pump (Watson Marlow 530U). The tank is filled step by step by first filling 5 cm bottom layer of well spread GB and then water was injected until the porous medium was fully saturated. This procedure should be carried out with great patience to avoid any air bubbles being captured inside the porous network. These bubbles, especially when they are connected to electrodes or thermocouples, can cause significant measurement errors. This operation is repeated step by step until reaching the top of the 2D tank and when all electrodes were buried by GB and water (almost 27.5 cm). A camera is placed on a fixed tripod in front of the 2D tank. Photos were taken during increasing CT levels and imbibition at regular time intervals of three and one minutes, respectively.

Previous experiments using imaging, time-domain reflectometry (TDR), and mass balance analysis showed a 2D flow behavior [START_REF] Colombano | Comparison of thermal and chemical enhanced recovery of dnapl in saturated porous media: 2D tank pumping experiments and twophase flow modelling[END_REF]. However, the character of the electric field can be in 3D, and therefore, our numerical model uses a 3D geometry (with a small thickness) to simulate electrical fields.

Experimental and measurement procedure Twenty-seven SIP measurements of an equatorial dipole-dipole (EDD) array were performed to resolve the variations of the electrical resistivity within the 2D tank. The conventional location of measurement for an EDD array is in the middle of the rectangles in Figure 1c. The exact place of measurements with the EDD array is illustrated in Figure 2 with blue circles and the number of measurements. Before and after each experiment, measurements were made with tap water to examine the performance and accuracy of electrodes and calculate the geometric factors (K). The geometric factor for water is easy to obtain because the resistivity of water is homogeneous and can be measured by a portable conductivity meter and the resistance by SIP lab. Then we calculated the K factor from resistivity over resistance. To validate the work, the measured geometric factors should be close to each other. According to our experience, the non-similarity of the geometrical factors was mainly related to non-polarizable electrodes or their surroundings. A blank test (with water only) was done for at least 24 hours to have enough data to increase measurement accuracy. The measurement readings may have serious errors, because of air trapping inside the polarizable electrodes. With a simple test with only water, these deficit electrodes can be replaced. In the next step, the experimental geometric factors (K) must be compared with ones that will be estimated by a numerical model for all 27 SIP measurements.

After creating a stable porous medium saturated with water, CT was injected from the bottom of the 2D tank (using the same points of entry where water was injected). For each step of increasing DNAPL content, measurements were made until the DNAPL has been reached to the middle of the 2D tank at 15 cm. This procedure was repeated three times to obtain a horizontal interface between CT and water in the middle of the tank. In the experimental process, obtaining a perfect flat boundary with CT was not easy due to fingering and dissimilar distribution in the different parts of the 2D tank. After stabilizing the experimental setup with the portion of half GB+water and the other half GB+CT, SIP measurements were made before heating, in the frequency range of 0.183 Hz to 20 kHz. Resistivity and temperature were measured at all levels in the 2D tank.

After heating the tank from the left side and diffusing the heat to the right side of the tank, our experimental protocol was done for drainage and imbibition of water where SIP measurements were conducted at each step. The heating time takes almost 12 hours to have a stable temperature in each thermocouple (measurement zone) of the tank. When the recorded temperature was stable, the pumping of CT from the bottom middle of the tank was started.

Due to the rapid imbibition, it was not possible to measure more than two frequencies during this dynamic stage of the experiments. Therefore, during the imbibition period, the SIP spectra were recorded in two frequencies (1.46 Hz and 2.92 Hz), close to the ones used commonly in the field scale. The resistivity calculations were, therefore, based on a single-frequency approach.

Iravani et al. (2020a) reported the α values (equation ( 2)) for the same liquids (water and CT) with and without GB in 1D cells and columns. From their comparison between the α value of water with other studies, normally the α value should be almost the same. We calculated this value again using our 2D tank experimental data to be sure it was the same in this new 2D configuration. In three levels of measurements, the top level is filled with

GB+water and the bottom level with GB+CT. Figure 3 is the linear regressions of the resistivity data with the temperature where the y-axis is the ratio of the resistivity at reference temperature (T 0 =20°C) and resistivity at T ( ) minus 1, and the x-axis is the temperature difference ( ). The slopes of the linear regression of the 9 measurements at the top with

GB+water and the bottom GB+CT are presented in Figure 3a and Figure 3b. The obtained regression is in accordance with the previous finding by Iravani et al. 2020a,b. Hayashi (2004) studied the linear relationship between electrical conductivity and temperature that was usually used in previous studies on water with different salinity and chemical substances. Regardless of the compositional and structural differences, the relationship between these two parameters of all different water samples was linear. He suggested being more cautious if temperature values are dissimilar to the reference temperature of 25 °C. The α values for water-saturated GB (0.052) and CT saturated GB (0.075) are the same as what we found in columns (Iravani et al., 2020a) but the coefficients of determination are less (GB+water 13.09% and GB+CT 3.28%) than columns for both.

According to Iravani et al. (2020b), the final equations for correcting resistivity for both multiphase systems are:

(6) (7)

After validating α values in two different geometries and despite the complexity of CT components, these two equations can be used in our 2D tank for temperature correction of electrical resistivity in the multiphase system of GB+CT and GB+water. The porous medium in the tank was subjected to heating with an immersion heater. The temperature fields (back view) in the tank before and after imbibition are shown in Figure 4. We used Surfer software to interpolate the data. Surfer interpolated the temperature and the corrected resistivity in the range of 110 to 100,000 Ω.m with the Kriging method. In this study, for all data fitted with Surfer, the same griding methods as Kriging were used with the experimental variograms fitted with linear component type and no trend removal. The temperature diffuses faster in the water at the top due to its higher thermal conductivity. The temperature close to the heater is almost 50°C and due to heat loss, it decreases across the tank to room temperature (20°C).

Experimental electrical resistivity value was measured and corrected with temperature for three testing EDD array levels. In the scatter plots of Figure 4, each point corresponds to a mean value of a square using an EDD array. Before imbibition, the first level of EDD at the bottom of the tank had almost the same resistivity. The decreasing resistivity in the middle part of the tank after imbibition suits the CT saturation decreasing and replacing DNAPL with water in this part. The points in the boundary of the tank were considered as same as their neighbor points for interpolation. Our study of temperature and saturation effects has been done separately for the condition that both factors want to be considered. After correcting the data for the temperature effect, the impact of saturation changes could be determined. The next step was pumping CT out from the middle bottom of the tank. Performing the imbibition step was to simulate the thermally enhanced pump and treat technology for soil remediation. The objective was to monitor how water and CT saturation vary inside the 2D tank using the electrical resistivity measurements.

Validation of the electrical resistivity measurements

Before starting the main drainage-imbibition experiment, to validate the performance of the electrodes, the first experiment was carried out by increasing the CT level in the tank at the normal room temperature (20°C). In the following, this experiment's data will be compared to the results obtained from a numerical model. After filling the half bottom of the tank with GB+CT, the resistivity was measured for 24 hours, then the temperature of the 2D tank was increased to 30 and 40°C by circulating hot water at a given temperature around the 2D tank to keep the temperature constant in all 2D tank parts. Compared to heating elements, heating with this method can result in temperature homogeneity throughout the whole tank.

For the main experiment, a laboratory protocol of a multiphase system was defined in the 2D tank that will be briefly described in the following. The 2D tank was filled with GB and two different liquids, water as a conductive fluid and CT with low conductivity. At the starting point, the tank was filled with GB and water layer by layer (Figure 5a). Figure 5b is the tank filled with GB and water, ready for injecting CT from the 5 points at the bottom of the tank (Figure 5c). It is necessary to slowly inject CT to have the relatively horizontal CT levels (Figure 5d). Despite all these considerations, in one case (Figure 5e) when a bump shape CT surface appears, slight injecting from the sides of the tank can balance the surface of CT in the tank.

We used an open-source modeling tool BERT (Boundless Electrical Resistivity Tomography) developed by [START_REF] Günther | Boundless Electrical Resistivity Tomography BERT 2-the user tutorial[END_REF] to model this simple drainage experiment. The possibility of modeling closed geometries as laboratory cells and tanks encouraged us to use this modeling code. The theory of developing this code based on the finite element method has been expanded from the previous study done by [START_REF] Rücker | Three-dimensional modelling and inversion of dc resistivity data incorporating topography-I[END_REF]. In the other studies, code was developed for arbitrary electrode shapes [START_REF] Rücker | The simulation of finite ERT electrodes using the complete electrode model[END_REF] and long electrodes (Ronczka, [START_REF] Günther | Boundless Electrical Resistivity Tomography BERT 2-the user tutorial[END_REF].

The 2D tank, by using three series of EDD arrays in different levels, was modeled by BERT. The fundamental physical theory for electrical resistivity used in BERT software is based on Ohm's law [START_REF] Günther | Three-dimensional modelling and inversion of DC resistivity data incorporating topography -II. Inversion[END_REF]. Assuming that the electric fields are conservative and the magnetic effects are negligible, the mathematical equations can be written as follow , ( 8)

where (A.m -2 ) is the current density, (V.m -1 ) is the electric field intensity, σ (S.m -1 ) is the electrical conductivity (the reciprocal of resistivity ρ), and (V) is the electric potential (primary variable). The initial condition was ρ(t=0)=30 Ohm.m that corresponds to the electrical resistivity of glass beads saturated by water. Neumann boundary condition (nocurrent flow boundary condition) was imposed on all 2D tank borders. The number of irregularly structured meshes was 1036942 elements.

In the other experiment in the 2D tank, increasing CT level is also simulated centimeter by centimeter from the bottom of the 2D tank up to 18.0 cm from the 2D tank bottom. Before these steps, we simulate a homogeneous tank filled with only GB and water.

After injecting CT, the resistivity values were recorded for these different levels, until reaching the top of the middle EDD level (17.5 cm). This data allows us to study how resistivity varies in the experiment and the model according to an increase in CT level.

Results and discussions

Correction of data with temperature

In remediating polluted sites with the thermal enhanced method, the effect of temperature on the electrical resistivity should be corrected before studying the effect of saturation change on geophysical parameters. At the beginning of the experiment, the mean temperature in the tank was 20°C. The temperature varied between 19.1 and 22.8°C. The resistivity was corrected with temperature by using the relation proposed by [START_REF] Dakhnov | Geophysical well logging: The application of geophysical methods; electrical well logging[END_REF] and [START_REF] Keller | Oxford. a case study in sichuan province, southwest china[END_REF] to have a geophysical parameter that is only controlled by saturation.

Electrical resistivity fields

The electrical resistivity variation was used to predict the water/CT saturation change.

The evolution of logarithmic apparent and inverted electrical resistivity (to facilitate the comparison before and after temperature correction) during the imbibition test is shown in Figure 6. The objective is to show how temperature changes can impact the calculation of resistivity. Iravani et al. (2020a) presented how the temperature affects the resistivity value in a multiphase sample. To study the effect of saturation change due to pumping, it is needed to deduct the effect of the heating. Figure 6a,b show the log of resistivity before temperature correction, and Figure 6c,d show the same results with temperature correction. This is the starting point of our study on the saturation effect on resistivity values in the 2D tank. Figure 6 is cut 10 cm from the left side of the tank to eliminate the influence of the metallic heater on the resistivity of the medium.

Validation of the setup configuration and electrical resistivity measurements using

BERT modeling

In modeling a physical phenomenon, the main constraint for a modeler is dealing with boundary conditions [START_REF] Stroud | Boundary representation modelling techniques[END_REF]. The first boundary problems are the Neumann boundary conditions [START_REF] Broggini | Immersive boundary conditions: Theory, implementation, and examples[END_REF]. We need a dipolar current injection to avoid overflow from the boundaries of the 2D tank [START_REF] Günther | Boundless Electrical Resistivity Tomography BERT 2-the user tutorial[END_REF]. Using a reference electrode inside the tank could solve this problem. This reference electrode should be placed slightly away from the current and potential electrodes inside the tank. The second problem that arises for

Neumann domains is the non-uniqueness of the solution of the partial differential equation due to a lack of calibration [START_REF] Günther | Boundless Electrical Resistivity Tomography BERT 2-the user tutorial[END_REF]. We used a calibration electrode to overcome this problem, which must be placed at the boundaries of the tank. The potential at this point (electrode) must be zero. The method used in the modeling was moving both electrodes (reference and calibration electrodes) in the tank and checking if the results obtained were consistent or not [START_REF] Günther | Boundless Electrical Resistivity Tomography BERT 2-the user tutorial[END_REF].

To validate our experimental measurements with numerical modeling, we completed two steps: In the first step, the ratio of geometrical factors K obtained from the experiment with water on the 2D tank and model were compared in Figure 7. In both the experiment and model, the K varies around 0.45 to 0.5 for the EDD arrays with ±20% error. We noticed that in points close to the metal grid on both sides of the tank (e.g., points 23, 24, and 25) the difference in K factors obtained from the model and experiment are more significant because, in the model, the effects of metal grids have not been considered.

In Figure 8, after increasing the level from 2 cm (entering the influence zone of EDD measurements) resistivity was increased both in the experiment and model. Since the influence zone of all measurements begins at 2.5 cm. In this figure, the model simulation results

properly coincide with the experimental data (Pearson's correlation coefficient, R=0.981). This validates the setup configuration and shows that the electrodes measure the resistivity with good precision. In the transaction part at a height of 3 and 4 cm, an overestimation was seen compared to the experimental resistivity data. This overestimation could be due to the effect of increasing CT and its direct effect on the resistivity value predicted by the model.

Phase variation during imbibition

During the imbibition period, both amplitude and phase shift of complex resistivity measured in the middle of the tank were affected by temperature and saturation changes.

Figure 9 shows how phase shift varied with time in the pumping influence zone in the middle of the tank for two different frequencies (f=1.46 and 2.92 Hz). Nine measuring zones (we call them points in the following) in the same vertical positions were chosen from three levels of the EDD array in the middle of the tank.

Points 12, 15, and 18 are at the top, consisting of the water-saturated zone (GB+water).

The results showed that the phase shift is close to zero in these three points. These results confirm previous findings on 1D cells and columns for the water-saturated cases (Iravani et al., 2020b). Points 11,14,17,and 10,13,16 are at the middle and bottom of the tank with an EDD array, respectively. During the imbibition (pumping process), CT was pumping from the middle bottom of the tank, and therefore, water replaced CT forming the cone of depression.

As shown in Figure 9, imbibition increased and then decreased phase shift at points 11, 14, and 17 (at the middle of the tank), but it was constant at points 10, 13, and 16 (at the bottom of the tank). This variation occurred sooner at point 14 than at point 11, which corresponds with our expectation because imbibition firstly affected the middle points, then the side points because of the cone of depression. Besides, the phase variations at points 14 and 17 happened at the same time. Drops in phase shift at the end of imbibition happened later for points 11 and 17 compared to 14 due to the delay of the imbibition effect from the middle of the tank reaching the sides. After imbibition (dominating water saturation), the phase shift dropped to a value close to zero. An increase in phase was observed in all 3 points in the middle of the tank, which can show the complex behavior of the spectra and phase due to the instantaneous change of water/CT saturation. However, we observed a weak response to imbibition tests for points 10, 13, and 16 situated at the bottom of the tank because they were not influenced by imbibition.

After dominating water in each block of the measurement, the recorded values became stable and close to zero. Comparing results of the phase change between Figure 9a and Figure 9b shows that small changes in frequency have a large effect on phase shift in layers containing DNAPL. The higher frequency, the higher the phase shift.

Predicting saturation in 2D tank

In a previous study, Iravani et al. (2020b) found a relation between saturation state and the amplitude and phase shift of the SIP spectra for CT+water, to have a calibration relationship for further uses in the 2D tank (current study) and even thereafter in field studies as below (Iravani et al., 2020b):

+ (9)
was obtained from the inverse relationship with the measured resistivity of the multiphase system in each zone. The only variable in equation ( 9) is water content , which can be present as water saturation after using this equation (S w = Ф). analysis to obtain the saturation. Therefore, we could not make a quantitative comparison of the saturation using these images.

Figure 10 shows resistivity measurements in different influence zones with interpolation that can give us the saturation field. In Figure 10c, we expected more saturation changes before and after the pumping at the left side of the tank compared to its right side. The heating from the left side has a direct effect on the viscosity of CT and consequently, the pumping rate, but the DNAPL/water interface is not horizontal and we could not quantify and make this assumption. As observed in Figure 6d and Figure 10c, the effect of imbibition seems to be well observed particularly in both sides compared to the middle.

The CT saturation was not exactly 100% even before starting imbibition because some percentage of water was captured in the pore spaces before injecting CT. It seems the Archie model does not have enough precision (the effect of surface conductivity should not be neglected) in very high resistivity (high CT saturation) and that monitoring in a small 2D tank does not allow us to capture the cone of depression with high precision for imbibition.

Effect of temperature and saturation on Cole-Cole parameters

In the non-isothermal condition before starting imbibition when the bottom half of the tank was filled with GB+CT, the complex resistivity was measured in the frequency domain (e.g., under electrical impedance conditions) before and after heating. The tank was heated from the left side by a heater and heat was dispersed from the left side of the tank to the right.

The most influenced part is close to the heater. There the temperature changed from 20 to 50°C, and that change affects two EDD measurements at the far-left side of the tank (No. 1 and 4). Figure 11 shows how the four Pelton parameters vary at 20 and 50°C. The experimental data were fitted to the Pelton model by using a Matlab code by [START_REF] Maineult | Upscaling of spectral induced polarization response using random tube networks[END_REF]. The fitting curves of the experimental data and model are presented in Figure 12.

We have already studied how these parameters vary with temperature in cells and columns (Iravani et al., 2020a). Studying these parameters in the 2D tank showed that  0 and log (τ) were decreasing with temperature increases. Raising the temperature from 20 to 50°C slightly increased c but both m and c were so close to 1 that it can be considered as a constant value. The fact that c is close to 1 means that the SIP response obeys a Debye model, and the fact that m is close to 1 means that the medium is highly chargeable because of poor dispersion of the spectra. These variations in results were the same for both measurements (1 and 4). Our findings on how the four parameters changed in the 2D tank confirm the previous findings on cells and columns with this difference: m was 0.89 and in the 2D tank for both measurements was almost 1.

Conclusions

The main goal of this study was to move one step closer to execute the usage of electrogeophysical methods to monitor the saturation changes in soils contaminated by a DNAPL (like CT) under non-isothermal conditions. In particular, applications of the SIP method in contaminated porous media were examined experimentally in the laboratory by using a 2D sample tank. Glass beads of 1 mm diameter were chosen because of similar physical properties to field soil. This allows to model as much as possible actual field conditions in the laboratory.

After preparing the non-isothermal condition by heating from the left side of the tank and correcting the temperature effect, we studied how the electrical complex resistivity varies during imbibition. Changes in the electrical complex resistivity magnitude were measured and the saturation field was obtained from the resistivity data. From the electrical resistivity measurements, we predicted the changes of the saturation field before and after imbibition accordingly. However, the proposed electrical resistivity monitoring did not produce the cone of depression during imbibition with high precision. This may be related to CT's very high resistivity and the small laboratory setup. The phase changes in 9 points in the middle of the tank in all three levels showed that the phase significantly decreases after increasing water saturation (imbibition period).

It is concluded that the temperature (5 to 7%) and saturation change due to replacing CT with water affects measured electrical complex resistivity. The effect of saturation change is much more considerable compared to the temperature. In the coupling of temperature and saturation variation after correcting resistivity with temperature (with a linear relation), the saturation effect will be the only variable in the real case with remediation surveys under the thermal-enhanced technique.

The effect of saturation change by injecting CT from the bottom of the tank was evaluated by BERT as a numerical modeling tool. The results show that the measurements with EDD arrays have been properly fitted by the BERT model. This validates our setup configuration and confirms that our electrodes functioned properly. The goal of our study was not to present comprehensive mechanistic modeling of the results but to present the experiment and the results that show the effect of heating and pumping on the SIP response, and its interest for field measurements. Developing a mechanistic model is a step beyond, the subject of future research.

The reported results showed that the proposed generalized Archie's law from the previous study could predict the saturation. However, during pumping at the bottom of the tank, the proposed coefficient did not perfectly reflect the saturation change. Finding the proper coefficients of generalized Archie's law can lead to precise saturation estimation in both high and low resistivity. These results at the bottom level during imbibition could be due to placing the current electrodes in the part that is saturated with CT. Our previous finding in 1D cells supports also the idea that in a case where the current electrodes are located in a saturated DNAPL part, SIP reports very high resistivity although the potential electrodes are still placed in the water-saturated part.

How Cole-Cole parameters change before and after heating were evaluated for 20 and 50°C before and after heating. The results of this study show that Cole-Cole parameter variations are in accordance with our previous findings.

We believe that the most outstanding part of our findings is the development of knowledge on how to monitor imbibition processes where both temperature and saturation change together in a 2D system. Validating our results in electro-geophysical field studies would advance their reliability in various environmental conditions. Our findings might be used for soils with similar physical properties.

FIGURE LEGENDS

Figure 1 Schematic layout of the a) front, b) top and, c) back of the 2D tank with the description of each element. In the front view of the tank, a double glazed transparent glass is used as a wall to better control the temperature and imaging purpose. Isothermal system includes a system of heating all around the tank using hot water. A stainless-steel heating element was inserted inside the left-side of the tank to produce the non-isothermal condition.

In sub-figure c) the red square showed an EDD measurement. There are 27 EDD measurements in three different levels. GB+CT for the bottom level of EDD array. The obtained α values for GB+water and GB+CT are in accordance with the results of the previous study on 1D columns (Iravani et al., 2020a).

Figure 4 Back view of the temperature field in the tank under non-isothermal conditions a) before, and b) after imbibition (the plus (+) signs are the positions of the thermocouples and red dot is the pumping point at the middle bottom of the tank). A significant temperature (54 and 18 °C) difference between the right and the left side of the tank is due to using a heater at the right side of the tank (in the back view). 1, 4, 7, 10, 13, 16, 19, 22, and 25; see Figure 2).

Figure

  Figure 10a,b, and Fig 10c,d show the saturation status and front image of the

Figure 2

 2 Figure2Equatorial dipole-dipole arrays from the front side of the tank. 10 current (metal) electrodes and 10 potential (non-polarized) electrodes in each level can record 9 EDD measurements. In three levels of measurements, 27 measurements can be recorded. Using other arrays (e.g., Wenner Beta) is also possible between two levels.

Figure 3

 3 Figure 3 Linear regression to find the α values for a) GB+water in the top level and b)GB+CT for the bottom level of EDD array. The obtained α values for GB+water and GB+CT are in accordance with the results of the previous study on 1D columns(Iravani et al., 2020a).

Figure 5

 5 Figure5Different steps of experimenting on the 2D tank. a) filling the tank with GB and water, b) tank packed with GB and fully saturated by water, c) injecting CT from 5 points at the bottom of the tank (for the places of injecting points, please see Figure1a), d) increasing step by step the level of the CT in the tank, e) possible bump shape of the CT surfacethe flat interface is ideal, f) almost half of the tank filled with GB+CT. This stage is the starting point of the imbibition experiment.

Figure 6

 6 Figure 6Log of the resistivity field in the tank a) before imbibition, b) after imbibition without temperature correction, and c) before and d) after imbibition with temperature correction, e) inverted log resistivity before temperature correction, and f) inverted log resistivity after temperature correction (the circle (•) and blue arrow signs are respectively the positions of the EDD measurements and drainage points).

Figure 7

 7 Figure 7 Experimental and model ratio of K factor for 27 measurements of electrical resistivity with EDD array. Point values closer to 1 show a better model prediction of the experimental data. Points far from the walls have better results because of the metallic grid used in the side walls.

Figure 8

 8 Figure 8Variation of resistivity as a function of CT height for EDD array at the bottom level of the tank. These are apparent resistivity obtained from the average values of the 9 measurements at the bottom level of the 2D tank (measurements No: 1, 4, 7, 10, 13, 16, 19, 22, and 25; see Figure2).

Figure 9

 9 Figure 9Variation of phase shift during imbibition period at the middle of the tank in 9 points (12, 15, and 18 at the top level, 11, 14, and 17 at the middle level, and 10, 13, and 16 at the bottom level) for the frequencies of a) 1.46 Hz and b) 2.92 Hz.

Figure 10 a

 10 Figure10a) The measured CT saturation field before imbibition, b) front image of the 2D tank before imbibition, c) the measured CT saturation field after imbibition, and d) front image of the 2D tank after imbibition. The black lines in "a" and "c" are the counters between two areas (GB+water and GB+CT).

Figure 11

 11 Figure 11 Effect of temperature on Cole-Cole parameters (20 and 50°C) close to the heater. As it was shown in Figure 5, just points close to the heater have recorded the temperature around 50 °C.

Figure 12

 12 Figure 12 Fitting curves of experimental data and model for a) point 1 at 20°C, b) point 1 at 50°C, c) point 4 at 20°C, and d) point 4 at 50°C. For more information about choosing the presented points, see Figure 11 description).

Acknowledgments

This study was performed as part of the Bioxyval project. The authors would like to thank ADEME (French Environment and Energy Management Agency) for co-funding the project under the "AMI" program. We gratefully acknowledge the financial support provided to the PIVOTS project by the "Région Centre -Val de Loire" and the European Regional Development Fund.

Data Availability Statements

The datasets generated and analyzed during the current study are available from the corresponding author on reasonable request.

Kemna, A., Binley, A., Cassiani, G., Niederleithinger, E., Revil, A., Slater, L., Williams, K.H., Orozco, A.F., Haegel, F.H., Hoerdt, A. and Kruschwitz, S. (2012) An overview of the spectral induced polarization method for near-surface applications. Near Surface Geophysics, 10(6), 453-468. Kremer, T., Schmutz, M., Leroy, P., Agrinier, P. and Maineult, A. (2016a)

Declaration of interests

☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:


Simon Benaïchouche
email: simon.benaichouche@imt-atlantique.fr

Guillaume Morel
email: guillaume.morel@imt-atlantique.fr

Franc ¸ois Rousseau
email: francois.rousseau@imt-atlantique.fr

Ronan Fablet
email: ronan.fablet@imt-atlantique.fr

Divergence-free continuous normalizing flows for uncertainty quantification

Uncertainty quantification in ill-posed inverse problems is a critical issue in a variety of scientific domains, including among others signal processing, imaging science, geoscience, remote sensing.... This has led to a variety of approaches, especially using Bayesian schemes such as Kalman methods, particle filtering schemes and variational Bayesian inference. Dealing with non-linear and non-Gaussian processes remain however a challenge, especially when considering high-dimensional systems. Recently, normalizing flows using deep neural networks have emerged as a very powerful tool to train generative models, which can sample realistic states, while making feasible the computation of the likelihood function. Their exploitation for uncertainty quantification in terms of differential Shannon entropy however requires the use Monte Carlo methods whose computational cost can be prohibitive, especially for high-dimensional systems such as time, space and space-time processes. Here, we introduce a new class of continuous normalizing flows with a divergence-free constraint for the underlying governing ordinary differential equations. This divergence-free constraint results in the preservation of the differential Shannon entropy through the trained flows. We demonstrate the relevance of the proposed framework to reach state-of-the-art performance for generative modeling tasks. We also illustrate applications to uncertainty quantification for the reconstruction of 1D and 2D states from partial observations. We discuss further our main contributions and applications to real-world case-studies.

Introduction

The reconstruction of hidden states from partial and/or noisy observations, referred to as inverse problems, are key challenges in numerous domains such as signal processing, computational imaging, medical imaging, remote sensing, geoscience,..... Model-driven schemes have long been the stateof-the-art methods, especially when dealing with physical states [START_REF] Evensen | Data Assimilation[END_REF]. They state the inverse problem as the minimization problem with two terms, a data-fidelity or observation term J and a re regularization term R:

X = arg min X J(X, Y ) + R(Y ) (1)
with Y the observation and X the state to be reconstructed.

The observation terms involves the definition of the forward observation model to relate X and Y . Regularization term states some prior on considered problem. For example when dealing with physics-driven dynamical systems, prior R is stated as R(.) := T k=1 ∥M(x k-1 ) -x k ∥ where M denotes a forecast operator associated to a physical equations [START_REF] Courtier | Variational Assimilation of Meteorological Observations With the Adjoint Vorticity Equation. Ii: Numerical Results[END_REF]. Recently, data-driven and deep learning methods have arised as relevant alternatives to address inverse problems. While deep learning approaches lead to the state-of-art-performance for inverse problems [START_REF] Dong | Learning a deep convolutional network for image super-resolution[END_REF], [START_REF] Minaee | Image segmentation using deep learning: A survey[END_REF] [START_REF] Xie | Image Denoising and Inpainting with Deep Neural Networks[END_REF] [START_REF] Yu | Generative image inpainting with contextual attention[END_REF], they also compete with model-driven schemes for the reconstruction and forecasting of geophysical systems [START_REF] Fablet | Learning Variational Data Assimilation Models and Solvers[END_REF][START_REF] Nonnenmacher | Deep Emulators for Differentiation, Forecasting, and Parametrization in Earth Science Simulators[END_REF].

These method are generally trained by gradient descent to minimize a given reconstruction cost.

θ * = arg min E(d(µ θ (Y ), X))

(2) it is known [START_REF] Bishop | Pattern recognition and machine learning[END_REF]) that when 2 is treated as a regression problem with least square estimation, a minimizer of 2 is given by the expectation of the conditional density p(X|Y ). The estimation of p(X|Y ) is of a key importance for in order to transpose many applications of a Data assimilation task in a stochastic perspective : One would like to take decision on the basis of a risk minimization for a wide range of applications instead of taking a decision relying only on the estimated mean µ θ . Furthermore, the knowledge of the probability density p(.|Y ) allows the computation of information theoretic related quantities such as the differential entropy : this quantities extends the classical Shannon entropy [START_REF] Shannon | The Mathematical Theory of Communication[END_REF] to continuous variable. For a random variable X with associated density function p, the differential entropy of X is defined by:

H(X) := E p [-log(p(x))] = - Ω p(x)log(p(x)).dx (3)
It can be understood as a measure of uncertainty : a random variable with high concentration of the measure would trend to have low entropy. One can give two examples for 1d variables : dirac distributions δ (no uncertainty) for which

H(δ) = -∞ or uniform distribution U on Ω = [a, b],
which corresponds to a case of maximal uncertainty verifies H(U ) = ln(b -a) → +∞ if |b -a| → +∞. and is of a key importance for applications such as sensor placement [START_REF] Yin | Entropy-based optimal sensor placement for model identification of periodic structures endowed with bolted joints[END_REF], However for an arbitrary distribution f , the computation of the differential Shannon entropy relies on Monte Carlo methods [START_REF] James | Monte Carlo theory and practice[END_REF] whose computational cost may be prohibitive for real-world processes.Thus derivation of approximation method of the differential entropy is an active area of researches [START_REF] Huber | On entropy approximation for Gaussian mixture random vectors[END_REF][START_REF] Hyvärinen | New approximations of differential entropy for independent component analysis and projection pursuit[END_REF][START_REF] Ao | Entropy estimation via normalizing flow[END_REF].

Here, we address uncertainty quantification according to the differential Shannon entropy in inverse problems, while keeping analytical an computation of the differential Shannon entropy to overcome the shortcomings of Monte-Carlo schemes for high-dimensional systems. We state the proposed framework as the learning constrained continuous normalizing flows [START_REF] Chen | Neural ordinary differential equations[END_REF]. Our main contributions are as follows:

• We introduce a novel class of divergence-free continuous normalizing flows. We show that the divergence-free constraint guarantees the preservation of informationtheoretic quantities such as the differential Shannon entropy. • We present a neural implementation of the proposed scheme with an explicit divergence-preserving parameterization of the neural architecture. • We asses the performance of the proposed scheme for generative modeling tasks with respect to state-of-the-art schemes and for application to uncertainty quantification in inverse problems.

This paper is organized as follows. Section 2 introduces the background and related work. We present for the proposed divergence-free continuous normalizing flows in Section 3. Section 4 introduces the resulting neural architecture and associated learning scheme. We report numerical experiments in Section 5 and further discuss our main contributions in Section 6. We include Proofs in Appendix.

Background and related work

Inverse problem and uncertainty quantification:

As stated in the introduction, we address uncertainty quantification issues in inverse problems. Let us assume that we are provided with a reconstruction scheme µ(Y ) for state X given partial observation Y . As mentioned above, we may consider both model-driven and learning-based schemes for the design of µ(•). Let ε(X) be the reconstruction error:

ε(X) := Y -µ(X) (4)
We aim at characterizing the density function p(.|y) associated to the reconstruction error ε p(.|y) and it's associated differential Shannon entropy. When considering Gaussian approximations for posterior p(x|y), one can derive analytically this quantity. It is however widely acknowledged that this approximation does not apply to most real-world processes which involve non-Gaussian statistics [START_REF] Evensen | Data Assimilation[END_REF]. State-of-the-art schemes also exploit Monte Carlo methods when one can sample posterior p(x|y) and compute log-likelihood log(p(ε(Y ))) The latter opens the floor to the exploitation of learning-based generative models such as normalizing flows [START_REF] Rezende | Variational inference with normalizing flows[END_REF].

Learning-based inference and normalizing flows:

From a probabilistic viewpoint, the maximum likelihood estimation (MLE) is the classic formulation for inference of probabilistic models. Given a parametric family of distributions according to parameter set Θ, we aim at identifying parameters θ * which maximize the likelihood of the reconstruction error ε(Y ) conditionnally to observation Y . This leads to the following minimization issue:

θ * = arg max θ∈Θ N i=1 log(p(x i -µ(y i )|y i , θ)) (5)
Deep learning has recently widened the class of parameteric probabilistic models one can consider in MLE framework through the introduction of normalizing flows. Normalizing flows [START_REF] Rezende | Variational inference with normalizing flows[END_REF][START_REF] Kobyzev | Normalizing flows: An introduction and review of current methods[END_REF]) make use of the formulas the change of variables to extend the MLE framework beyond classic parametric probabilistic families. Given a latent random variable u of density p, the density p of the image v = f Φ (u) is given by:

p(v) = 1 |Jf Φ (u)| p(u) (6)
where Jf Φ (u) denotes the Jacobian matrix of f Φ evaluated at point u. Given a parameterization of the likelihood of the latent variable, usually as a Gaussian distribution, normalizing flows then come to design mapping operator f Φ so that one can compute the determinant of its Jacobian. Among normalizing flow approaches, continuous normalizing flows (CNF) [START_REF] Chen | Neural ordinary differential equations[END_REF] states mapping operator F Φ as a flow operator governed by an ordinary differential equation:

∂u(t) ∂t = f Φ (u(t)) (7)
Under this parameterization, the resulting mapping for likelihood values involve the divergence of the flow:

∂log(p(u(t)) ∂t = -T r ∂f Φ (u(t)) ∂u(t) = div (f Φ ) (8)
This equation relates to a Lagrangian point of view of Liouville equations which describe the infinitesimal evolution of the log-density by a flow governed by the ODE ∂u ∂t = f Φ (u(t)). This framework allows to fit complex distributions without computing the determinant of ( 6) which is O(N 3 ). The price to pay is the numerical integration of the ODE.

Equation ( 8) also stresses that divergence-free flows satisfying div(f Φ ) = 0 would even simplify the computation of likelihoods to the evaluation of the likelihoods in the latent space. We explore this avenue in this work and show it also leads to an analytical derivation of the differential Shannon entropy for image variable v. While the parameterization of divergent-free flows is relatively straightforward for 2D and 3D processes using Helmoltz decomposition of velocity fields, we also address higher-dimensional processes. We may also point out that divergent-free flows also relate to some extent to Halmitonian flows. This class of flows is governed by some underlying energy functions. Recent exploration of trainable Halmitonian flows [START_REF] Greydanus | Hamiltonian neural networks[END_REF][START_REF] Sanchez-Gonzalez | Hamiltonian graph networks with ode integrators[END_REF] for the data-driven discovery of dynamical systems emphasize the relevance of constrained neural ODE schemes to regularize learning problems.

Divergence-free continuous normalizing flow

Here, we focus on a particular class of continuous normalizing flow by enforcing a divergence-free constraint on the CNF. This constraint relates to an incompressibility or volume-preserving constraint on the vector fields defined by the CNF. We show that the proposed divergencefree CNF allows us to derive an analytical computation of information-theoretic quantities such as the differential Shannon entropy. Formally, let us consider the following CNF:

     u(0) ∼ N (µ, σ 2 ) ∂u(t) ∂t = f Φ (u(t)) s.t. div [f Φ (u(t))] = 0 (9)
Early litterature on normalizing flow already investigated the use of volume preserving flow [START_REF] Dinh | Nice: Nonlinear independent components estimation[END_REF] with additive transformation whose triangular matrix is triangular making the computation of its determinant trivial. the proposed approach allows the computation of more complex flows.

Parameterization of divergence-free flows The construction of the particular class of divergence-free normalizing flow presented here rely on the ability to represent divergent-free functions. Let us introduce the following parameterization:

f Φ = A.∇Ψ (10) where A is an antisymmetric linear operator satisfying A T = -A and Ψ a learnable scalar function. Overall, parameters Φ combine the parameterization of matrix A and of scalar function Ψ. As detailed below, functions of the form are divergent free. For a more intuitive physical interpretation of the presented approach, we adopt a graph representation of the function described above.

Graph representation of divergence-free flows: For any point of R n , any governing ODE x ′ (t) = f (x(t)) is described by a graph (V, E) and a potential function Ψ. For any edge (x i , x j ) ∈ E, we associate the stream function defined as

S ij (x) = [0, ..., - ∂Ψ(x) ∂x j , 0, .., ∂Ψ(x) ∂x i , 0, ...0] (11)
whose i -th component equals -∂Ψ(x) ∂xj and j -th component equals ∂Ψ(x) ∂xi . For each (x i , x j ) ∈ E we then associate a weight w ij , the dynamical model f is stated as : The graph (V, E) describes the coupling between data features which in turn relates to ODE x ′ (t) = f (x(t)). The following lemma state that f can be split onto independent dynamic corresponding to components of (V, E) which are divergent-free:

f (x) := (i,j)∈E w ij S ij (x) (12)
Lemma 1. Let assume that the graph (V, E) is composed of N components (V n , E n ). Then the equation x ′ (t) = f (x(t)) can be split into N coupled non-autonomous equations x ′ n (t) = F n (x(t)) with F n verifying : T r(∇ X k F k (x)) = 0 (13)
As a corollary, we retrieve the expected divergent-free feature for the considered parameterization. Corollary 2. The function f described in (Eq.11) is divergent-free i.e : div(f ) = 0 (14)

From this property, we can derive the volume-preserving property of the considered CNF: Lemma 3. Let f be a C 1 (R n , R n ) function. Then, for any t ∈ R + , the semi-group Γ(t) of solutions of the ODE :

x ′ (t) = f Φ (x(t)) is volume preserving i.e : |JΓ(t)(x)| = 1 (15)
for any x ∈ R n where JΓ denotes the jacobian of flow Γ.

Robustness to mode collapse:

Training continuous normalizing flows under the maximum likelihood framework may be subject to generate probability density function with singularities. Let think about a dynamical system f for which each training point x i is attractive. This would result in the negative likelihood function going to +∞. From a dynamical system perspective, a point x i is attractive if f (x i ) = 0 and if the eigenvalues of the linearized system have negative real-part eigenvalues. As stated by the following theorem, the considered CNF has no attractive points, which in turn should prevent from mode collapse issues.

Theorem 4 (Robustness to mode collapse). Let x ∈ R n be a zero of the divergent-free function described above i.e f Φ(x) = 0. Then the eigenvalues of Df (x) have zero real parts.

Derivation of information-theoretic quantities

We further exploit the key features of the proposed CNF to derive information-theoretic quantities. More specifically, the flow defined by (Eq. 9) preserves the differential Shannon entropy as a direct by-product of the divergence-free constraint.

Lemma 5. Let be q a C 1 (R n , R) scalar function, and Q the scalar quantity defined as :

Q(p) = R n q(p(x))dx (16) If Φ ∈ Dif f (R n ) denotes a diffeomorphism of R n which satisfies |JΦ(x)| = 1.
Then, the push-forward density p := p#Φ verifies

Q(p) = Q(p) (17)
In particular, p and p have same differential Shannon entropy.

As the divergence-free constraint implies the volumepreserving property, the later applies to the proposed CNF to derive information-theoretic quantities. More specifically, the flow defined by (Eq.9). A direct consequence of this lemma is that the differential Shannon entropy of any push-forward density image according to (Eq.9) is solely determined by the entropy of initial distribution N (µ, σ 2 ) given by H = 1 2 ln{(2πe) d det(σ 2 )}. As such, it avoids computationally-expensive computation required by classical Monte-Carlo methods which integrates the highdimensional ODE for each sample.

Figure 2: Samples generated from a diagonal gaussian after 2 days of training : The divergence-free flow was limited to a 8-neighbours graph structure, it fail to generate realistic samples, however we discuss in appendix of way to improve the expressiveness of the presented approach.

Learning framework

Given an arbitrary algorithms µ(.) which delivers a reconstruction of state X from partial observations Y , we address the characterization of the reconstruction error using a neural implementation of the proposed divergence-free CNF. The proposed implementation relies on CNN architectures satisfying the divergent-free constraints presented above. As stated in (Eq.11), our model relies on learning a CNN parameterization for potential Ψ and a linear antisymmetric transformation. We first describe the latter. We then introduce the resulting end-to-end neural architectures and the associated training setting.

The antisymmetric transformation layer: In order to apply the proposed framework to high-dimensional ndimensional space, we propose to implement a translationinvariant antisymmetric transformation in a CNN-fashion. This can simply be done using convolutional transform and applying a antisymmetric operator to the channel dimension.

as described in Figure 3.

Figure 3: For 1-channel image, the antisymmetric transform layer consist in 3-step.

Step 1 : a fixed convolution layer is applied to an input, which are stacked onto 4 channels.

Step 2 : a antisymmetric operator is applied over the channel.

Step 3 : the inverse of operation performed in step 1.

Neural architecture : overall, the considered neural architecture for the modeling and characterization of the posterior p(X|y) is sketched in Fig. 4. The considered neural architecture involves two main components:

• A first neural network Θ which takes as inputs the observation Y and the reconstructed state µ(X) and outputs the initial condition of the flow. It combines a prediction of the the covariance of the initial Gaussian distribution N (µ, σ 2 ) (we consider a diagonal coavariance parameterization); z 0 : an augmented component used by the neural integrator to better advect the initial gaussian density. We may recall that the estimation of the mean µ(x) is given by an arbitrary denoising algorithm µ. • The second neural network consists in a Runge-Kutta integration scheme ofof the considered divergent-free neural flow f Φ applied to initial condition u(0) = (ε(X), z 0 ). For the numerical integration, we make use of the open source library torchdiffeq provided by the authors of [START_REF] Chen | Neural ordinary differential equations[END_REF].

Training setting :

Given the considered neural architecture, we train the proposed divergence-free CNF using a MLE criterion (Eq.??) as considered in [START_REF] Chen | Neural ordinary differential equations[END_REF][START_REF] Grathwohl | Ffjord: Free-form continuous dynamics for scalable reversible generative models[END_REF]:

L(θ) = ln(|σ|) + (z i -µ(x)) t σ -1 (z i -µ(x)) (18)
with z i := Φ(-t)(y i , Θ(x)) the output of the flow integrated backward in time. We also perform training of a gaussian model with full covariance for benchmarking purpose.

All experiments were run using pytorch. The training procedure is performed by gradient descent using Adam algorithm over 150 epoch. As mentioned above, the neural solver for ODE is the pytorch implementation of [START_REF] Chen | Neural ordinary differential equations[END_REF].

Numerical experiments

This section reports numerical experiments for the proposed approach. We first assess divergence-free CNF for generative modelling tasks before considering applications to inverse problems and uncertainty quantification.

Benchmarking for generative modeling tasks : We perform density estimation tasks using the benchmarking framework proposed in [START_REF] Papamakarios | Masked autoregressive flow for density estimation[END_REF]. The reported quantitative comparison involves state-of-the-art normalizing flow approaches, namely the masked autoregressive flows (MAF) [START_REF] Papamakarios | Masked autoregressive flow for density estimation[END_REF], the masked autoencoder for distribution estimation (MADE) from [START_REF] Germain | Made: Masked autoencoder for distribution estimation[END_REF] and the continuous normalizing flow. All schemes are evaluated in terms of log-likelihood of the test set. We report state-of-the-art performance, better than MADE [START_REF] Germain | Made: Masked autoencoder for distribution estimation[END_REF], with a clear improvement compared with the Gaussian baseline. Given that the proposed scheme can be regarded as a constrained version of [START_REF] Grathwohl | Ffjord: Free-form continuous dynamics for scalable reversible generative models[END_REF], we do not expect to outperform [START_REF] Grathwohl | Ffjord: Free-form continuous dynamics for scalable reversible generative models[END_REF] in general. Sampled examples illustrated in Fig. 2 for MNIST dataset may in this respect point out that we may increase the complexity of the our implementation to sample more realistic examples. We may however stress that the considered parameterization was primarily considered here for an application to unvertainty quantification as illustrated below.

Reconstruction of Lorenz-63 dynamics

We illustrate the relevance of the proposed framework for uncertainty quantification in inverse problems with dynamical systems. As toy model, we consider Lorenz-63 system which involve chaotic dynamics:

x ′ 1 (t) = σ(x 2 -x 1 ) x ′ 2 (t) = x 1 (ρ -x 3 ) -x 2 x ′ 3 (t) = x 1 x 2 -βx 3 (19)
We asses the ability to estimate the conditional probability p(x T |y 1:T ) of the state x T given all previous noisy observations Y 1:T -1 . We state the observational model as

H(x 1 , x 2 , x 3 ) = x 1 + x 2 .
Given a time series of noisy observations, a LSTM (µ) is trained to reconstruct the true state while an other LSTM (Θ) is trained to estimate the parameters of the flow : In (de Bezenac Emmanuel 2020) authors propose the same framework for filtering, the estimation of the conditional distribution was performed using realNVP [START_REF] Dinh | Density estimation using real nvp[END_REF]. We report results in Table 2 where we compare the learned model with a Gaussian baseline. They clearly point out that the proposed scheme outperforms a Gaussian approximation for the posterior, which is the classically considered including in learning-based variational Bayesian scheme [START_REF] Kingma | Auto-encoding variational bayes[END_REF]. Table 2: Reconstruction performance in terms of conditional log-likelihood over the test set for Lorenz-63 case-study: we compare the proposed approach to a Gaussian baseline. We refer the reader to the main text for the experimental setting. Figure 6: Differential entropy of the estimated distribution p(ε|y 1 ) over the test dataset :As expected, the entropy of the estimated density is higher over the middle region where it is hard to predict in which part of the attractor the true state lies due to the chaotic nature of the governing equations.

MNIST case-study : We consider an inpainting application for MNIST dataset and aims at assessing how the observation patterns impact the reconstruction performance. We set two observation configurations with the same number of observed pixels set to d = 81: the first configuration involves regularly-sampled pixels in both horizontal and vertical directions, the second configuration samples pixels regularlyspaced along the horizontal axis and randomly along the ver-tical axis. We also add a gaussian white noise to the data to ensure that the reconstruction error does not lie in an area with zero Lebesgue measure. We trained two neural network µ 1 , µ 2 on the image space to perform the reconstruction from the gappy observations. For each reconstruction network, we train a divergence-free CNF to model the reconstruction posterior and evaluate the associated Shannon entropy. We reduce train and test sets to digits from 7 to 9 in order to reduce training time and power consumption of the experience. We perform training over 600 epoch using Adam algorithm and stop each models after 36 hours of training, no improvement over the training set were observed. We report in Fig. 8 examples of the random sapling patterns, of mean reconstruction and samples from the inferred posterior for the randomly-sampled observations. We report in Fig. 9 the distributions of reconstruction error of the inferred Shannon entropy for the two observation configurations. The random sampling configuration leads on average to smaller reconstruction errors, our model infers more similar distribution of reconstruction uncertainty, with even a greater uncertainty for the iregular sampling configuration. This is in line with compressive sampling results. While random sampling strategies improve almost surely the reconstruction performance compared with a regularly-spaced sampling, it may also lead to reconstruction outliers, which in turn results in heavier-tail for the posterior distribution.

We interpret the results reported in Fig. 9 as an illustration of the proposed scheme to capture differences in the non-Gaussian tails of the posterior through the inferred Shannon entropy. how volume-preserving transformations and the underlying graph structure of the dynamical system affects the learning process.

Figure 1 :

 1 Figure 1: An example of graph (V,E) with two components : nodes represents the data features, edges their connection through the existence of a function S ij described below. In terms of physical interpretation the existence of an edge ij authorize the learned model to advect the measure in the plan of R 2 corresponding to features x i , x j

Figure 4 :

 4 Figure4: Given an arbitrary denoising algorithms µ, the proposed neural network architecture aim to learn the density function p associated to the reconstruction error ε : first a CNN neural network Θ (green) extract both features θ and log-variance ln(σ) of an initial gaussian distribution, then a divergence-free flow parametrized by θ advect the gaussian distribution using an ODE solver (red) to match the reconstruction error distribution ε in the sense of maximumlikelihood estimation.

  Figure7further reveals the non-Gaussian features of the posterior when the state is close to the bifurcation zone from one lobe to another one of Lorenz-63's attractor. This Figure plots samples of the inferred posterior p(X T |Y 1:T ) for different values of T for a true state close to the bifurcation zone. As expected, the more observations, i.e. T being larger, the lower the spread of the samples, which indicates a lower uncertainty. Figure5further illustrates this global pattern through the comparison of the inferred Shannon entropy averaged over the test set compared with the mean square error of the reconstruction. We may recall that we cannot derive analytically the true posterior, nor true Shannon entropy for Lorenz-63 dynamics. We further illustrate the relevance of the proposed derivation of the Shannon entropy of the posterior in Figure6. The inferred entropy alongside the attractor of the Lorenz 63 nicely emphasizes that the reconstruction uncertainty is larger when getting closer to the bifurcation zone.

Figure 5 :

 5 Figure 5: Evolution of the mean differential Shannon entropy associated to the inferered posterior p(X T |Y 1::T ) with respect to time horizon T for Lorenz-63 case-study (red). For comparison purposes, we consider the evolution of the mean square error (MSE) of the reconstruction of the true state (blue).

Figure 7 :

 7 Figure 7: Left : true and inferred Lorenz-63 states and states, right : sampling of 1000 points of the inferred posterior p(x T |y 1:T ) for different value of T. It clearly illustrates the ability of the proposed framework to address non-Gaussian features.

Figure 8 :

 8 Figure 8: left : masked input y, center : estimated reconstruction µ(y), right : a data sample s generated using the learned distribution ε by the additive relation s = µ(y) + ε

Table 1 :

 1 Caption : performances in terms of likelihood over the test set (higher is better).

	method	POWER GAS HEPMASS MNIST
	Gaussian	-7.74	-3.58	-27.93	-1366
	MADE	-3.08	3.56	-20.98	-1380
	FFJORD	-0.46	8.59	-14.92	NA *
	MAF	0.14	9.07	-17.7	-1300
	VPCNF (ours)	-1.96	6.32	-20.13	-1176

* : authors provides results in terms of bits per dim

Appendix

proof of lemma 1 Suppose the graph (V, E) is composed of N components (V n , E n ) : then the equation x ′ (t) = f (x(t)) can be split onto N coupled non-autonomous equations x ′ n (t) = F n (x(t)) with F n verifying :

T r(∇ X k F k (x)) = 0 (20)

Proof. Without loss a generality, we suppose

The last line follow from the Cauchy Lemma

Proof of lemma 4:

Proof. if we denote as z = Φ(x), the change of variable formula states that p(z), the image density is given by p(z) = p(x)

It follows from dz = |JΦ(x)|dx and |JΦ(x)| = 1 that :

Thus, the differential Shanon entropy H(p) of an image density associed to a measure-preserving flow is invariant. It is entierly determined by the initial density p and ca be computed analytically if H(p) does.

Proof of theorem 4:

which is antisymetric as the matricial product between the antisymmetric matrix A and the symmetric hessian matrix ∇ 2 Ψ(x). Thus its eigenvalues have zero real part.

limitations of the proposed approach

Given a divergent free flow F and a familly P 0 of initial probability density, It may arise that the presented framework here still lack of expressiveness in order to fit an arbitrary density probability f . For example, the set of 1dimensional divergent free flow contains only the translations. Here we discuss about two ways to overcome this issue :

Graph structure and data representation : All experience performed in this paper (except for the Lorenz-96) make use of the same 8-neighbours graph-structure as described in figure 10. Because it remains fixed during the training, this may limit the performances of the proposed framework for density estimation tasks. In future works we may investigate the optimization of this graph and consider different representation of the input. In [START_REF] Voleti | Multi-Resolution Continuous Normalizing Flows[END_REF] authors provides an interesting multi-scale representation using volume preserving transformations. Augmented state Given a set of examples x i ∈ R n , we can construct an a set of augmented variable xi := [x i , ε i ] with ε i ∼ N (0, σI) and apply the above mentioned framework This may help to match an arbitrary distribution f with the marginal density of the learned probability density. The price to pay is that the Shanon entropy of the marginal density of f is no more given by the initial density, we can only compute a lower bound given the inequality :

Extend the set P 0 In this work we only considered diagonal gaussian distributions as set P 0 of initial density function, the action of volume preserving on this set results in a non-universal approximator. However, it extends any mathematical results on estimation of the differential Shanon entropy of any initial set P.