O Lengliné

J E Elkhoury

G Daniel

J Schmittbuhl

R Toussaint

J.-P Ampuero

M Bouchon

Interplay of seismic and aseismic deformations during earthquake swarms: an experimental approach

Observations of earthquake swarms and slow propagating ruptures on related faults suggest a close relation between the two phenomena. Earthquakes are the signature of fast unstable ruptures initiated on localized asperities while slow aseismic deformations are experienced on large stable segments of the fault plane. The spatial proximity and the temporal coincidence of both fault mechanical responses highlight the variability of fault rheology. However, the mechanism relating earthquakes and aseismic processes is still elusive due to the difficulty of imaging these phenomena of large spatiotemporal variability at depth. Here we present laboratory experiments that explore, in great detail, the deformation processes of heterogeneous interfaces in the brittle-creep regime. We track the evolution of an interfacial crack over 7 orders of magnitude in time and 5 orders of magnitude in space using optical and acoustic sensors. We explore the response of the system to slow transient loads and show that slow deformation episodes are systematically accompanied by acoustic emissions due to local fracture energy disorder. Features of acoustic emission activities and deformation rates distributions of our experimental system are similar to those in natural faults. On the basis of an activation energy model, we link our results to the Rate and State friction model and suggest an active role of local creep deformation in driving the

Introduction

Numerous observations of a correlation between seismic activity enhancements and slow slip transients in the Earth crust have been reported in various tectonic locations [START_REF] Linde | A slow earthquake sequence on the San Andreas fault[END_REF][START_REF] Crescentini | Constraints on slow earthquake dynamics from a swarm in central Italy[END_REF][START_REF] Lohman | Earthquake swarms driven by aseismic creep in the Salton Trough, California[END_REF][START_REF] Segall | Earthquakes triggered by silent slip events on K īlauea volcano, Hawaii[END_REF][START_REF] Liu | Seismicity variations associated with aseismic transients in Guerrero, Mexico, 1995[END_REF] and in geothermal areas [START_REF] Bourouis | Evidence for coupled seismic and aseismic fault slip during water injection in the geothermal site of Soultz (France), and implications for seismogenic transients[END_REF][START_REF] Takada | Aseismic slip during the 1996 earthquake swarm in and around the Onikobe geothermal area, NE Japan[END_REF]. Also, observations of postseismic slip and aftershocks following large earthquakes are manifestations of transient deformation coupled with abundant earthquake activity [START_REF] Perfettini | Postseismic relaxation driven by brittle creep: A possible mechanism to reconcile geodetic measurements and the decay rate of aftershocks, application to the Chi-Chi earthquake, Taiwan[END_REF]. A central question is whether or not these two different mechanical responses of faults, i.e. seismic and aseismic slip, occur on closely located zones. If they are, one expects a strong interaction between the two processes.

In some cases, the seismic signal concurrent with slow slip events is characterized as tectonic tremors or low frequency earthquakes, as observed in subduction zones [START_REF] Rogers | Episodic Tremor and Slip on the Cascadia Subduction Zone: The Chatter of Silent Slip[END_REF][START_REF] Obara | Episodic slow slip events accompanied by non-volcanic tremors in southwest Japan subduction zone[END_REF][START_REF] Ito | Slow Earthquakes Coincident with Episodic Tremors and Slow Slip Events[END_REF] or in transform tectonic settings like the San Andreas Fault (SAF) [START_REF] Vidale | A survey of 71 earthquake bursts across southern California: Exploring the role of pore fluid pressure fluctuations and aseismic slip as drivers[END_REF][START_REF] Nadeau | Nonvolcanic tremors deep beneath the San Andreas fault[END_REF][START_REF] Shelly | Migrating tremors illuminate complex deformation beneath the seismogenic San Andreas fault[END_REF]. The aseismic nucleation phase of mainshocks has also been related to foreshocks, as in the case of the 1999 Izmit earthquake [START_REF] Bouchon | Extended Nucleation of the 1999 M w 7.6 Izmit Earthquake[END_REF] or the 2011 Tohoku-Oki earthquake [START_REF] Miyazaki | Seismic and aseismic fault slip before and during the 2011 off the Pacific coast of Tohoku earthquake[END_REF][START_REF] Ando | Possibility of Mw 9.0 mainshock triggered by diffusional propagation of after-slip from mw 7.3 foreshock[END_REF][START_REF] Kato | Propagation of slow slip leading up to the 2011 mw 9.0 Tohoku-Oki earthquake[END_REF]. Seismic events associated with aseismic slip are generally located on fault planes suggesting that they represent a dynamic shear instability on the sliding interface (La [START_REF] Rocca | Cascadia Tremor Located Near Plate Interface Constrained by S Minus P Wave Times[END_REF][START_REF] Shelly | Precise location of San Andreas Fault tremors near Cholame, California using seismometer clusters: Slip on the deep extension of the fault?[END_REF]. Geodetic inversions suggest that the aseismic motion occurs on the same fault plane as the seismic events [START_REF] Lohman | Earthquake swarms driven by aseismic creep in the Salton Trough, California[END_REF]. However the deformation signal recorded at the surface by GPS and/or InSAR instruments only provides a macroscopic view of the deformation process at depth. Details of the aseismic slip distribution are often lacking due to the limited resolution, unless the aseismic slip occurs at the Earth surface [START_REF] Doubre | Fluid-controlled faulting process in the asal rift, djibouti, from 8 yr of radar interferometry observations[END_REF].

The close spatial and temporal occurrence of both seismic and aseismic slip suggests a causal relation between the two phenomena. Nonetheless, the causal mechanism is not straightforward as earthquakes can both trigger and be triggered by slow slip events [START_REF] Du | Triggered aseismic fault slip from nearby earthquakes, static or dynamic effect?[END_REF][START_REF] Perfettini | Postseismic relaxation driven by brittle creep: A possible mechanism to reconcile geodetic measurements and the decay rate of aftershocks, application to the Chi-Chi earthquake, Taiwan[END_REF]. It might also be that both aseismic and seismic slip are manifestations of a common deformation process. In this case, the seismic signal can be seen as the signature of patches on the fault plane deforming dynamically whereas the geodetic observation represents an integrated signal over all the deforming sites. Therefore, the deformation on a fault plane takes place over a wide range of speeds. This is supported by observations of heterogeneous postseismic slip on the SAF inferred from the analysis of repeating earthquake sequences [START_REF] Lengliné | Inferring the coseismic and postseismic stress changes caused by the 2004 M w = 6 Parkfield earthquake from variations of recurrence times of microearthquakes[END_REF]. The dynamic events represent the high velocity tail of the slip speed distribution while its average produces the observed geodetic signal.

The contribution of earthquakes to the total amount of slip released during transient episodes is generally small [START_REF] Lohman | Earthquake swarms driven by aseismic creep in the Salton Trough, California[END_REF]. Despite being located on a common interface, accurate spatial location of the seismic activity relative to the aseismic slip is difficult to obtain. The heterogeneity of mechanical and physical properties of the interface might control the partition between seismic and aseismic slip, the interface being envisioned as a collection of brittle patches (asperities) embedded in an otherwise creeping region which accommodates slow slip [START_REF] Linde | A slow earthquake sequence on the San Andreas fault[END_REF][START_REF] Lohman | Earthquake swarms driven by aseismic creep in the Salton Trough, California[END_REF][START_REF] Wech | Seismic and geodetic constraints on Cascadia slow slip[END_REF][START_REF] Perfettini | Seismic and aseismic slip on the Central Peru megathrust[END_REF]. The concentration of asperities, or brittle patches, governs the relative importance of dynamic failures in the deformation process. It is also readily possible that temperature produces a broad scale effect on the slip partition over the interface and thus constrains the transition between brittle and ductile rheology.

The physics of both seismic and aseismic processes is not straightforward as it involves a complex problem on a spatially heterogeneous medium with a large number of degrees of freedom and short and long range interactions. Several numerical models have tried to reproduce the evolution of such systems but computations are extremely demanding and time consuming and limited to only large scale heterogeneities [START_REF] Kaneko | Towards inferring earthquake patterns from geodetic observations of interseismic coupling[END_REF], or to the quasi-dynamic approximation [START_REF] Hillers | Statistical properties of seismicity of fault zones at different evolutionary stages[END_REF][START_REF] Ariyoshi | Migration process of very low-frequency events based on a chain-reaction model and its application to the detection of preseismic slip for megathrust earthquakes[END_REF].

Analog laboratory experiments overcome these difficulties as space-time integrations are performed by the deforming system without any model requirement. Here, we investigate the slow (creep) and fast (acoustic) deformations produced by the propagation of an interfacial brittle-creep crack. Our experimental setup is a much simpler configuration than that of a complex fault zone system but the fundamental processes of interest, the interplay of slow deformation and brittle fracture on the same heterogeneous interface, remain similar. Our setup allows the simultaneous monitoring of both acoustic activity and the detailed geometrical evolution of the fracture. We show that seismic and aseismic events co-exist in the system and their activity rates are highly correlated. This suggests that seismic activity during swarm episodes is driven by the elastic loading of asperities by local aseismic deformations.

Experimental Setup

Sample Preparation

To prepare each sample, we use two transparent PMMA (poly methyl methacrylate) plates of dimensions 20 × 10 × 1.0 cm and 23 × 2.8 × 0.5 cm (Figure 1). First, we sand-blast one surface of the narrow plate with glass beads of diameter φ ∈ [180 -300] µm. We clean the blasted plate to remove any electrically attached glass beads. Then we assemble the two plates in a stiff metallic loading cell with the blasted surface facing one of the surfaces of the thick plate. Finally, we impose a homogeneous normal load on the assembled plates and heat the loaded sample to 190 • C for 45 minutes to anneal the plates. The thermal annealing produces a cohesive interface, weaker than the bulk, along which the sample will break under load. The sand-blasting introduces random roughness to the plate surface that controls the heterogeneity of local strength along the interface. It also induces microstructures on the plate surface that make the surface opaque. The newly formed block, after annealing, recovers its transparency since the contrast of the refraction index along the interface disappears (see [START_REF] Grob | Quake catalogs from an optical monitoring of an interfacial crack propagation[END_REF]; Lengliné et al. (2011) for details). Interestingly PMMA exhibits a brittle behavior at short time scales and is semi-brittle or even plastic at longer times.

Macroscopically this long time scale regime is described by a ductile rheology. PMMA shows a time-temperature equivalence desirable for addressing either high temperature processes or very long term evolution [START_REF] Ward | An introduction to the mechanical properties of solid polymers[END_REF]. This richness of the PMMA rheology enables the observation of a brittle-creep rupture regime. It provides an attractive analogy for the study of numerous time-dependent mechanisms in the Earth crust as those originating at the brittle-ductile transition.

Mechanical loading

Once the sample is ready, we clamp the widest PMMA plate to a stiff aluminum frame. A stepping motor applies the loading over the top side of the narrow plate in a direction normal to the plate interface (Figure 1).

We measure the vertical displacement of the loading point with a linear variable differential transformer (LVDT) with a resolution of 1.3 µm. The vertical displacement imposed on the narrower plate induces stable mode I propagation of a planar fracture along the prescribed weak interface. We impose a variation of the loading speed, u(t), to simulate a variation of the external driving force. We performed 14 experiments where we applied various forms of temporal transients, namely a step increase, a bump increase and a sinusoidal of the loading speed imposed by the stepping motor. All these transient loads produced broad scale variation of the front speed (e.g.

Figure 2, see also Supplementary Material) and they could be interpreted as an analog for far field stress perturbations caused, for example, by pore fluid pressure, magmatic intrusion or mantle flow.

Optical events

We monitor the fracture front propagation using a fast camera (Cam-Record 600) with up to 1000 fps or a slow speed camera (Nikon D700) with up to 5 fps to follow the progression of the front position over longer time scales (Figure 1). Optical images of the interfacial rupture show dark and bright regions corresponding to open crack and unbroken parts of the sample respectively. Image processing determines the transition between dark and bright areas that defines the fracture front. We first compute the difference between each image and the first image of the experiment. The image difference highlights the fracture front while removing permanent artifacts. Then, grayscale images are transformed into black and white images according to a gray level threshold separating bright and dark regions. Then, we calculate the gradient in the direction of front propagation to highlight the transition zone. We finally extract connected pixels from the gradient images that correspond to the front position, a(x, t). The front propagates along the y axis with the origin defined at the loading point and is positive in the direction of crack propagation. The x axis is perpendicular to y and defines the coordinate of a point along the front and ā(t) is the mean position of the front at time t (see [START_REF] Måløy | Local waiting time fluctuations along a randomly pinned crack front[END_REF] and [START_REF] Grob | Quake catalogs from an optical monitoring of an interfacial crack propagation[END_REF] for details).

We compute the local speed of the crack as the pixel size divided by the time spent by the front in each pixel (Figure 3). We extract the pixels with the highest speeds from the local movements of the front. The broad scale variation of the front velocity could be compared to slow events recorded by geodetic instruments in a tectonic context. The patches of higher than average deformation speed could be considered as analogous to local creeping episodes or slow slip events. They represent slow deformation episodes that locally exceed the macroscopic deformation rate. We follow the procedure detailed by [START_REF] Grob | Quake catalogs from an optical monitoring of an interfacial crack propagation[END_REF] to define such events. We define an "optical event" (OE) as a as a cluster of connected pixels which velocity exceeds a prescribed speed threshold, v th . Here we adopt v th = 10 × v , where v is the mean velocity of the crack during an experiment. Although these optical events have higher velocities than most of the deforming sites, they remain slow deformation episodes compared to dynamic events (the average speed v ≃ 500µms -1 , far from the Rayleigh wave speed, V r ≃ 1.7 10 9 µms -1 ).

These optical events are characterized by a Gutenberg-Richter power law with a slope similar to that of tectonic earthquakes [START_REF] Grob | Quake catalogs from an optical monitoring of an interfacial crack propagation[END_REF]. This Gutenberg-Richter relation might also be related to the magnitude-frequency scaling inferred for slow-slip events [START_REF] Wech | An earthquake-like magnitude-frequency distribution of slow slip in northern Cascadia[END_REF].

Acoustic events

The crack propagation produces acoustic activity that we monitor with a 32 elements linear array of piezo-electric sensors. Sensors are located on a line parallel to the plate axis and sensors separation is 3 mm (see Fig. 1).

The typical distance between the closest acoustic sensor and the border of the plate is 1 cm. The sensors peak frequency response is ∼ 500 kHz and all channels are continuously recorded at 5 MHz. For each experiment, we extract the acoustic signal recorded in the two closest acoustic sensors to the fracture front line. We manually trigger the acquisition of the acoustic data, which is synchronized with the camera time sequence. The recordings of the acoustic signals last, in our experiments, for a maximum of 50 s.

We use a short-term-average to long-term-average ratio (STA/LTA) to detect acoustic events (AE) in the recorded acoustic signal [START_REF] Earle | Characterization of global seismograms using an automatic-picking algorithm[END_REF]. This procedure is similar to that applied to earthquake data. In our case, we use shorter time windows tailored to our short signal durations and high frequency acquisition rate. The STA window length is 20µs, the LTA window length is 100µs and the threshold for setting an STA/LTA detection is 2.0. An event must be detected simultaneously on the two closest channels to be considered in our analysis. We show an example of the acoustic signal recorded during one experiment in Figure 4. The typical duration of the recorded events is on the order of 100µs. The performance of the detection algorithm is illustrated in Figure 5. As expected, the detection of events is associated with high amplitudes of the recorded signal on the two closest channels. AEs result from the fracturing of geometrical asperities over the interface. We estimate the typical spatial scale of these failures to be less than several microns (The largest AE is generated during an optical event that has maximum dimensions of ∼ 600 µm by 30 µm). We observe nonuniformly distributed clustering of acoustic and optical activity (Figure 5).

The temporal coincidence of AEs and OEs is important but not systematic (Figure 5). It is attributed to AEs that are too small to be optically detected, OEs that are genuinely aseismic, and clusters of AEs that are lumped into a single large OE.

Results

Distribution of inter-event time

Our results suggest the presence of clustering of the acoustic activity (Figure 5). We quantify this clustering by calculating the distribution of inter-event times of successive AEs for all 14 experiments. We only considered AE when the loading rate was nearly constant in order to avoid mixing populations recorded during different loading rates. Inter-event times are normalized by the average AE rate of each experiment (the average AE rate for the different experiments, is on the order of 10 2 events/sec). The probability density functions (pdf) for all the experiments are well approximated by a gamma function (Figure 6). The displayed gamma distribution is obtained from averaging the parameters of each individual fit. Figure 6 also includes the inter-event time distribution for the ISC worldwide catalog for earthquakes with magnitude M>5.5 in the period 1975-2004. Temporal clustering has been well documented for earthquakes and is suggested to be a result of interactions among earthquakes [START_REF] Corral | Long-term clustering, scaling, and universality in the temporal occurrence of earthquakes[END_REF][START_REF] Molchan | Interevent Time Distribution in Seismicity: A Theoretical Approach[END_REF][START_REF] Hainzl | Estimating Background Activity Based on Interevent-Time Distribution[END_REF][START_REF] Saichev | Theory of earthquake recurrence times[END_REF]. The normalization of the interevent time is obtained similarly by the average seismic event rate (on the order of 1 event/day). The good agreement between the gamma distribution and the pdf obtained from our experimental data is similar to that for earthquake data [START_REF] Corral | Long-term clustering, scaling, and universality in the temporal occurrence of earthquakes[END_REF] or AE recorded during rock fracture experiments [START_REF] Davidsen | Scaling and universality in rock fracture[END_REF]. It suggests that the mechanism responsible for the seismic activity in fault systems or rock fractures shares strong similarities with our experiment, which can be seen as an analogous system for seismicity interactions. We interpret the decay at short time scale as the emergence of an Omori-Utsu law describing interactions among events while the distribution evolves to an exponential distribution representative of a Poisson process at longer time scales [START_REF] Saichev | Theory of earthquake recurrence times[END_REF].

Power-law distribution of local slip-rate

We calculate the pdf of the local velocities following [START_REF] Måløy | Local waiting time fluctuations along a randomly pinned crack front[END_REF] for six of the experiments where we used the fast video camera at high sampling rate. We only analyze, for these experiments, time windows with an almost constant large scale velocity. The local velocity is normalized by the average velocity of the crack, v . Front velocities exhibit significant fluctuations at small scales (Figure 3). For v > v the pdf shows a power law decay with exponent 2.55 consistent with previous works [START_REF] Måløy | Local waiting time fluctuations along a randomly pinned crack front[END_REF]Lengliné et al., 2011). This power law behavior of the deformation rate at small scale, in our experiment, has to be compared to the observed behavior in natural fault systems. For instance, in southern California, fault slip rates are found to obey such a power law scaling where the fault system is characterized by slow slip rates [START_REF] Meade | Power-law distribution of fault slip-rates in southern California[END_REF]. This comparison suggests that a second relation between our experiments and natural fault behavior can be proposed and that the scaling behavior is a general feature of slowly deforming media in the presence of heterogeneities.

Lateral migration of the deformation velocity

Migration of the fracture velocity along the crack front direction is analogous to the dip-parallel slip propagation revealed by tremors in subduction zone [START_REF] Shelly | Complex evolution of transient slip derived from precise tremor locations in western Shikoku, Japan[END_REF][START_REF] Ghosh | Rapid, continuous streaking of tremor in Cascadia[END_REF]. In most experiments, high loading rates or low recording acquisition rates do not allow us to distinguish migration patterns confidently as interactions among sites are numerous and occur over short time intervals. We find indications of along-front migration with tractable velocities for long enough times during experiments where loading rates are the lowest (Figure 8). High velocity clusters appear in streaks that suggest the progression of the rupture in the direction perpendicular to the crack front propagation. The migration velocity of these high velocity streaks is on the order of 2 cm/s (Figure 8). During this time interval, the crack front propagated at a typical velocity of 400 µm/s. We note that this migration speed is higher than the crack front speed. Taking 10-100 km/h as a typical range of dip-parallel tremor migration [START_REF] Shelly | Complex evolution of transient slip derived from precise tremor locations in western Shikoku, Japan[END_REF][START_REF] Ghosh | Rapid, continuous streaking of tremor in Cascadia[END_REF], and 10 km/day as an estimate of along-strike speed of the slow slip event (e.g. [START_REF] Bartlow | Space-time correlation of slip and tremor during the 2009 Cascadia slow slip event[END_REF], the ratio of these two velocities is in the range 25-250. The velocity ratio obtained in our experiment is on the order of 50, well within the range obtained in subduction zones.

An interplay of seismic and aseismic local deformations

The evolution of the acoustic emission (AE) rate, the spatial average velocity of the crack front and the optical event (OE) rate are similar (Figure 7). AE, OE rates and the average front velocities are computed for intervals of 0.2 s. A clear temporal relation between these three quantities at this large scale can be inferred: the AE and OE rates closely follow the macroscopic variation of crack front speed. The brittle fracture propagation generating acoustic emission and creep deformation co-exist in the same zone and are active at the same time.

We calculate the affine relation connecting the average front speed and the AE rate for all experiments. The residuals of the linear relation are small and symmetrically distributed suggesting a linear relation between AE rate and crack velocity. This is attested by the correlation coefficient, ρ, computed for each experiment between these two variables. The correlation coefficient is defined as ρ = σ xy /σ x σ y , with σ xy the covariance of the AE rate with the front speed and σ x and σ y the standard deviations of the AE rate and the front speed. The correlation coefficient, ρ for a 0.2 s window and for our 14 experiments ranges from ρ = 0.75 to 0.97 and with a mean value of 0.87. We now focus on smaller time scales relevant to acoustic activity and slow movement of the crack front. The crack front velocity shows important fluctuations at small scales, although being smooth and rather continuous when observed at large scale. These fluctuations result from the heterogeneous nature of the interface and the elastic interactions along the crack front line.

We investigate the temporal relation between OE and AE recorded during the 6 experiments which show the highest acoustic activity and the best resolved crack advance. The cross-correlation function, C(∆t), between the rate of AE, r AE (t), and the rate of optical events, r OE (t) is

C(∆t) = 1 √ C AE C OE t f t 0 r OE (t + ∆t) × r AE (t), (1)
where C AE and C OE are the auto-correlations values at zero lag time of the AE and the OE rates, respectively. For both AE and OE, rates are computed as the number of events detected in time intervals of 5 ms from time t 0 to time t f (the mean is removed from the time-series). The maximum correlation between the two types of signals is at zero lag time (Figure 9). It shows that the acoustic activity occurs mostly in a short time span around the aseismic deformation. The quasi-symmetric shape of the correlation function also suggests that optical events both precede and follow acoustic activity during the course of an experiment highlighting the close interplay between these two modes of deformation. The slight asymmetrical shape of C(∆t) (Figure 9) also suggests that a higher rate of OE occurs following AE than preceding it, similar to postseismic slip observed after large earthquakes (e.g. [START_REF] Hsu | Frictional Afterslip Following the 2005 Nias-Simeulue Earthquake, Sumatra[END_REF].

Discussion

In our experiments, the average fracture speed obeys the Arrhenius law at the macroscopic or global scale (Lengliné et al., 2011). Such an evolution has been inferred in numerous experimental systems reproducing slow deformations of rocks and other materials [START_REF] Atkinson | Subcritical crack growth in geological materials[END_REF]. In our experiment, this expression relates the crack front speed propagation, V , to an activation energy mechanism and can be expressed as

V = V 0 exp -Q k B T . ( 2
)
V 0 is the product of a frequency of attempt to break molecular bonds and the jump distance, k B is the Boltzmann constant and T is the temperature.

The activation energy of the process is

Q Q = Q 0 -Gα 2 , ( 3
)
where Q 0 is the activation energy in the absence of applied stress, G is the energy release rate and α is the typical length scale of the atomistic fracturing process. The stress dependence in Eq. ( 3) is similarly often written in terms of the stress intensity factor, K, at the crack tip. It follows that by controlling the imposed loading rate on the system we modify the value of the energy release rate (equivalently stress intensity factor) at the crack tip. In this sense, simulated transient deformation processes mimic the increase of the deformation speed recorded during swarm episodes [START_REF] Lohman | Earthquake swarms driven by aseismic creep in the Salton Trough, California[END_REF] or slow slip events in subduction zones [START_REF] Rogers | Episodic Tremor and Slip on the Cascadia Subduction Zone: The Chatter of Silent Slip[END_REF].

The formulation of the crack speed in Eq. ( 2) has been proposed to be the physical basis of the rate and state friction model widely considered in fault mechanics (e.g. [START_REF] Nakatani | Conceptual and physical clarification of rate and state friction: Frictional sliding as a thermally activated rheology[END_REF][START_REF] Rice | Rate and state dependent friction and the stability of sliding between elastically deformable solids[END_REF]. If we consider a volumetric process, the corresponding expression to Eq. ( 3) is

Q = Q 0 -τ c α 3 (4)
where the stress dependence is written in terms of the average shear stress, τ c along asperity contacts. α 3 is the activation volume of the process. Following [START_REF] Rice | Rate and state dependent friction and the stability of sliding between elastically deformable solids[END_REF], we call σ c the average normal stress on contact zones, defined as: τ c /σ c = τ /σ, with τ and σ being the large scale stresses. It follows from Eq. ( 2)

V = V 0 exp -Q 0 + τ c α 3 k B T (5)
leading to:

τ = σ Q 0 α 3 σ c + k B T σ c α 3 ln V V 0 ( 6
)
where V is interpreted in terms of slip rate on the fault plane. Eq. 6 has a similar form to the empirical relation of the rate and state friction model with a = k B T /σ c α 3 and the state dependent variable represented by Q 0 /α 3 σ c [START_REF] Rice | Rate and state dependent friction and the stability of sliding between elastically deformable solids[END_REF]. It follows that the progressive evolution of the front position in our system can be viewed as a similar mechanical problem to the evolution of slip on a fault plane with the same constitutive equation. The quasi-static evolution of an anti-plane shear rupture is governed by similar formulations to mode I fractures [START_REF] Gao | Shear Stress Intensity Factors for a Planar Crack With Slightly Curved Front[END_REF][START_REF] Schmittbuhl | Slow Crack Propagation and Slip Correlations[END_REF]. Thus, we consider our mode I fracture problem analogous to shear ruptures. Therefore our experimental setup is relevant for addressing the mechanics of simplified fault models as studied numerically by [START_REF] Gao | Penetration of a quasi-statistically slipping crack into a seismogenic zone of heterogeneous fracture resistance[END_REF]; [START_REF] Perfettini | Shear and normal load perturbations on a two-dimensional continuous fault: 1. Static triggering[END_REF]. Provided this analogy, the crack front advance a(x, t) is analogous to the slip along fault δ(x, t). However each point of the interface that breaks is completely unloaded due to the geometry of our experimental setting. Therefore, we have no repeating failure of the same patch contrary to natural faults and AE stress drops are total. However, earthquake stress drops are usually a small fraction of the total stress (e.g. [START_REF] Kanamori | Mechanics of earthquakes[END_REF]. Nevertheless, in both systems, AE and earthquake still represent a dynamic failure mode, which makes them comparable. Following the approach presented by [START_REF] Dieterich | A constitutive law for rate of earthquake production and its application to earthquake clustering[END_REF], we propose to relate the AE rate to the stress history of our system. We assume an interface composed of a population of sources that generate AE. We hypothesize that the duration of nucleation of these sources is longer than the duration of stress variations imposed to the system. Hence, the seismicity rate, R, varies exponentially with the shear stress change ∆τ imposed on the system [START_REF] Dieterich | A constitutive law for rate of earthquake production and its application to earthquake clustering[END_REF][START_REF] Beeler | Why earthquakes correlate weakly with the solid Earth tides: Effects of periodic stress on the rate and probability of earthquake occurrence[END_REF][START_REF] Cochran | Earth Tides Can Trigger Shallow Thrust Fault Earthquakes[END_REF][START_REF] Helmstetter | Afterslip and aftershocks in the rateand-state friction law[END_REF],

R = r exp ∆τ aσ (7)
where r is a reference seismicity rate. The analogy mentioned above between shear stress τ and the energy release rate G, suggests that in our experiment:

R = r exp α 2 ∆G k B T , (8)
Notice that Eq. ( 8) is the same expression of that for the speed of the crack as a function of G obtained by combining Eq. ( 2) and ( 3)

: V = V * 0 exp (α 2 G/k B T ) with V * 0 = V 0 exp (-Q 0 /k B T ).
Therefore we have a direct linear relation between AE rate R and speed of the crack V : R ∝ V , (Figure 7). At the system size, this relation provides a link between acoustic activity and periods of high deformation. At smaller scales, the same mechanism is applicable as the acoustic activity is influenced by local important deformations episodes (OE) (Figure 9). Such local creeping episodes are similar to the small localized creep events observed on the San-Andreas fault by [START_REF] Scholz | Seismic and Aseismic Slip on the San Andreas Fault[END_REF].

We emphasize that the deformation recorded at the surface of the Earth during transient episodes should be viewed as the macroscopic integration of sequences of creep events. On large spatial scales, the creep events formed the observed large scale behavior and correspond to the macroscopic observation of the transient slip speed variation. Such a view is similar to the concept developed by [START_REF] Scholz | Seismic and Aseismic Slip on the San Andreas Fault[END_REF] where the deformation is modeled as the integration of discrete sites subject to time-dependent strength. The acoustic deformation being recorded only represents the part of the deformation corresponding to zones that fail dynamically. This model is fully compatible with our observations. We also remark that sites do not fail independently as attested by the local close temporal activity between creep events and acoustic events. As proposed by [START_REF] Scholz | Seismic and Aseismic Slip on the San Andreas Fault[END_REF], it is readily possible that creep events recorded at our experimental scale are also formed of small acoustic events such that down-scaling might exist up to the microscopical level associated with the discrete breakage of single molecular bonds.

Conclusion

We analyze the coupled evolution of acoustic activity and slow deformation during the propagation of brittle-creep fractures in a heterogeneous medium. Our unique experimental setup addresses the relation between seismic and aseismic slip along natural faults. It provides an original characterization of slow deformation processes which are difficult to capture on faults at depth. Numerous statistical features of the deformation observed in natural systems are reproduced by our experiments like the gamma distribution of interevent times and the power law distribution of slip rates. We show that the acoustic activity (dynamic events) is part of the deformation process and it occurs over a widely distributed range of speeds including slow slip. Also, small scale observations of the deformation reveal coexistence between creep and acoustic events in the same mechanical system. The small scale complex dynamics lead to a macroscopic integrated signal of the deformation that shows a smooth and continuous deformation speed that correlates with the rate of acoustic events. (∆t) between the rate of acoustic events (AE) and the rate of optical events (OE). AE rate and OE rate are computed as the number of events during intervals of 5 ms and mean is remove from the time-series. The crosscorrelation function corresponds to an averaged function computed over 6 experiments. We observe that the maximum of the correlation function is found at zero time lag. We also notice the increase of the correlation function around the peak, suggesting that OE are clustered in time for some duration before and after an AE. The inset figure shows the normalized autocorrelation functions for the AE (black curve) and the OE (gray curve). Both functions show some increase around zero time lag supporting the interplay between these two modes of deformation.

> We build an original experiment designated to study the relation between slow and dynamic deformations in the brittle creep regime. > We observe numerous acoustic events in relations with local creeping episodes. > In relation with faulting processes, we propose that earthquake swarms are driven by local slow slips.

*Highlights

Figure 1 :

 1 Figure1: Left: Side view of the experimental setup. The bottom plate is separated from the upper one using a loading force applied by a rod connected to a stepping motor. The upper PMMA plate is attached to a stiff aluminum frame (short dashed lines). The load causes a deflection u of the bottom plate and the propagation of an interfacial crack. The crack front is located at distance ā from the loading point. The front advance is monitored by a high or slow speed camera set in vertical position, perpendicular to the crack plane. Right: detailed bottom view of the sample and the loading axis. The dark gray zone corresponds to the cracked zone. The multiple black areas show the acoustic sensors of the linear array.

Figure 2 :

 2 Figure 2: Evolution of the loading force (red), and the loading displacement (blue) during an experiment. The crack is supposed to start moving at t = 45s as evidenced by the force peak. The black box on the top figure represents a zoom displayed on the bottom figure. A transient variation of speed is superimposed to a constant loading displacement rate between t = 65 s to 85 s (the black line corresponds to the average front position). Camera and acoustic time windows are displayed respectively as dark gray and light gray shaded areas.

Figure 3 :

 3 Figure3: A: Map of the local speeds of the front. Image scale is given by the length of the vector showing the front propagation direction which is 1.3 mm long. The front propagates from bottom to top. We observe small scale fluctuations of the crack front speed. Black dots represent hypocenters of optical events obtained after thresholding the velocity map and are displayed at the centroid of the corresponding high velocity cluster. B: Zoom on a subzone of the interface located in the black rectangle in A. C: Probability density functions (pdfs) of the local velocities computed for 6 experiments. We observe a power law decay of the pdfs for v > v with an exponent ν = -2.55 compatible with[START_REF] Måløy | Local waiting time fluctuations along a randomly pinned crack front[END_REF].

Figure 4 :Figure 5 :

 45 Figure 4: Example of recorded signal of an acoustic event (AE). We observed a modification of the frequency content associated with the arrival of the AE wave train. The duration of the AE on the displayed channel is around 100µs. Sampling rate is 5 MHz

Figure 6 :

 6 Figure6: Distribution of normalized inter-event time for all experiments. Each experiment is represented by a different colour. The best gamma distribution fitting all the experiments, is represented by a black curve: p(τ ) = Cτ γ-1 e -τ /β , where C = 0.44, β = 1.8 and γ = 0.54. Black filled diamonds show the inter-event time distribution computed from the ISC worldwide catalog for earthquakes with magnitude M > 5.5.

Figure 7 :

 7 Figure7: Evolution of the rate of AE (purple diamonds), of the average crack front speed (blue squares) and of the optical event (OE) rate (red circles) as a function of time. Rate are computed for interval of 0.2 s for a time period encompassing the loading transient shown in Figure2. We observe that the rates of AE and OE are well correlated with the variation of the crack front speed at this broad time scale.

Figure 8 :

 8 Figure 8: Top: Crack front speed as a function of time and position along front during one experiment. High velocities appear as streaks that extend along the crack front direction. Bottom left: zoom that corresponds to the space time domain delimited by the black rectangle in the upper figure. White dashed lines show the migration during two high velocity episodes. Slope of these lines gives an estimated migration speed of the order of 2 cm/s.

Figure 9 :

 9 Figure9: Cross-correlation function C(∆t) between the rate of acoustic events (AE) and the rate of optical events (OE). AE rate and OE rate are computed as the number of events during intervals of 5 ms and mean is remove from the time-series. The crosscorrelation function corresponds to an averaged function computed over 6 experiments. We observe that the maximum of the correlation function is found at zero time lag. We also notice the increase of the correlation function around the peak, suggesting that OE are clustered in time for some duration before and after an AE. The inset figure shows the normalized autocorrelation functions for the AE (black curve) and the OE (gray curve). Both functions show some increase around zero time lag supporting the interplay between these two modes of deformation.

  Figure 1aClick here to download Figure: fig1a.eps

  Figure 2 Click here to download Figure: fig2.eps

  Figure 4 Click here to download Figure: fig4.eps

  Figure 5 Click here to download Figure: fig5.eps

  Figure 7 Click here to download Figure: fig7.eps

  Figure 8 Click here to download Figure: fig7bis.eps

  Figure 9 Click here to download Figure: fig8.eps

Acknowledgments

We thank K. J. Måløy, K. T. Tallakstad, S. Santucci, M. Grob, F. H. Cornet, J.P. Avouac, D. R. Shelly, M. Aktar and H. Karabulut for fruitful discussions. We also thank the editor and two anonymous reviewers for their suggestions and A. Steyer for technical support. We acknowledge the support of ANR grant SUPNAF and of NSF grant EAR-1015698.


Khalifa Nsir

Gerhard Schäfer
email: schafer@unistra.fr

Raphaël Di

Chiara Roupert

Olivier Razakarisoa

Renaud Toussaint

Laboratory experiments on DNAPL gravity fingering in water-saturated porous media

Keywords: gravity-driven instabilities, immiscible displacement, porous medium, DNAPL, optical fibers, residual saturation

Introduction

Characterization and quantification of the migration of immiscible liquids in groundwater are topics that have received considerable attention in recent years. Recent studies of two-phase immiscible flows have been motivated by the need to treat spills and leaks of socalled dense nonaqueous-phase liquids (DNAPLs), such as trichloroethylene (TCE), which can severely impact the quality of subsurface water supplies [START_REF] Kueper | An overview of immiscible fingering in porous media[END_REF][START_REF] Birovljev | Gravity invasion percolation in 2 dimensions-experiment and simulation[END_REF][START_REF] Fayers | On the importance of gravity and three-phase flow in gas displacement processes[END_REF][START_REF] Bettahar | Surfactant Enhanced In Situ Remediation of LNAPL Contaminated Aquifers: Large Scale Studies on a Controlled Experimental Site[END_REF][START_REF] Benremita | Quantification du transfert de trichloroéthylène en milieu poreux à partir d'un panache de vapeurs vers la nappe d'eau souterraine[END_REF][START_REF] Bohy | Caractérisation de zones sources de DNAPL à l'aide de traceurs bisolubles : mise en évidence d'une cinétique de partage[END_REF]. Unstable immiscible displacement can form viscous, capillary, and/or gravity fingers. These fingers propagate rapidly, causing early breakthrough relative to stable displacement [START_REF] Brailovsky | Fingering Instability in Water-Oil Displacement[END_REF].

The current physical understanding of fingering processes has been developed primarily through linear stability analysis [START_REF] Yao | Stability analysis of the unsaturated water flow equation 2. Experimental verification[END_REF]. While this theory provides information about the displacement process under unstable conditions and the practical growth rate of fingers for small amplitude perturbations, it neither predicts nor explains many significant features observed in laboratory experiments. The assumptions inherent to the linear analysis are somewhat limiting with respect to actual field conditions. Fingering processes were initially studied in 2D experiments, such as empty Hele-Shaw cells, as in the fundamental work of [START_REF] Saffman | The penetration of a fluid into a porous medium or Hele-Shaw cell containing a more viscous liquid[END_REF]. These cells are composed of two closely spaced parallel plates, and the motion of fluids in the system is mathematically analogous to that of twodimensional flow in a porous medium. Later on, porous media was added between the plates of Hele-Shaw cells for better reproduction of the random character of porous geometry.

Observations from laboratory experiments on instabilities have allowed researchers to both test the predictions of the linear stability theory for fingers and to discover a number of unexpected features of fingering flow fields, in which instabilities arise from the presence of gravity (Frett et al., 1997), viscosity [START_REF] Méheust | Interface scaling in a twodimensional porous medium under combined viscous, gravity and capillary effects[END_REF][START_REF] Toussaint | Influence of porescale disorder on viscous fingering during drainage[END_REF]Løvoll et al., 2010), both gravity and viscosity (Méheust et al., 2002, Nicholl and[START_REF] Nicholl | Experimental observations of gravity-driven fingering during immiscible displacements within a rough-walled fracture[END_REF] or capillarity when the displacement is very slow (Lenormond et al 1989, Cottin et al;2010).

Laboratory studies have also been conducted on soil columns to characterize the behavior of DNAPL displacement in saturated porous media and to better understand the complexities of DNAPL migration in environmental media [START_REF] Khataniar | The effect of reservoir heterogeneity on the performance of unstable displacements[END_REF][START_REF] Riaz | Influence of Relative Permeability on the Stability Characteristics of Immiscible Flow in Porous Media[END_REF][START_REF] Tullis | Wetting front instabilities: a three-dimensional experimental investigation[END_REF]. Schwille (1988) conducted large column experiments (1 m high, 0.4 m diameter) with perchloroethylene (PCE) in a low-permeability saturated sand (K=2×10 -4 ms -1 ). He found that the PCE front was not uniform or even bulbshaped. Instead, the PCE penetrated as tendrils or fingers in the interior of the sand body. Schwille (1988) also investigated PCE penetration into water-saturated glass beads with diameters ranging from 0.49 to 0.7 mm. The DNAPL was observed to penetrate into the saturated zone in a finger-like stream. Anderson et al. (1992) described the common occurrence of DNAPLs in the saturated zone of aquifers in the form of fingered plumes and pools that cause erratic contaminant distributions. [START_REF] Glass | Physics of gravity fingering of immiscible fluids within porous media: An overview of current understanding and selected complicating factors[END_REF] reported that a complicated pattern of interacting fingers resulted from their stratified soil experiments; and this limited the number of isolated fingers for which velocity and diameter data could be collected. To explore the wetting phase saturation levels in fingers, [START_REF] Geiger | Infiltration in homogeneous sands and a mechanistic model of unstable flow[END_REF] and [START_REF] Dicarlo | Experimental measurements of saturated overshoot on infiltration[END_REF] conducted one-dimensional fingering experiments using small cylindrical columns in which they examined the variation in DNAPL saturation along the finger length. Their results suggested that the propagating fingertip is near full saturation, and the saturation level decreases in the direction of the finger tail. [START_REF] Glass | Gravity destabilized non-wetting phase invasion in macroheterogneous porous media: Experimental observations of invasion dynamics and scale analysis[END_REF] used a high-resolution light transmission method to investigate gravity-destabilized non-wetting phase invasion (CO 2 and trichloroethylene) in a macro-heterogeneous porous medium. They observed a series of fingers and pools, behind the growing front, as a sequence of gravitystabilized and gravity-destabilized displacements within the heterogeneous medium.

To describe the two-phase liquid-liquid displacement in a porous medium, various dimensionless parameters such as the capillary number N c (ratio of viscous to capillary forces), the viscosity ratio M (ratio between the displaced fluid viscosity and the displacing fluid viscosity), and the Bond number N b (ratio of gravity to capillary forces) are commonly used in the literature. These numbers are defined as follows:

where v is the Darcy velocity [LT -1 ], μ w [ML -1 T -1 ] is the dynamic viscosity of the wetting fluid, σ ow [MT -2 ] is the interfacial tension between the nonwetting (oil) and wetting fluid (water) and θ [-] is the contact angle at the solid-water-DNAPL interface;

where μ 1 [ML -1 T -1 ] is the dynamic viscosity of the displacing fluid and μ 2 [ML -1 T -1 ] is the dynamic viscosity of the displaced fluid, and

where Δρ [ML -3 ] is the difference between the density of displacing fluid and the density of displaced fluid, g [LT -2 ] is gravitational acceleration, a [L] is the mean grain diameter of the porous medium, and  [-] is the angle formed by the flow vector and gravity vector. When

    cos v N ow w c (1) 2 1 M    (2) θ cos ) cos( 2 ow b ga N        (3)
gravity forces are involved, the Bond number can be either positive or negative, depending on the value of . The common convention is that N b is considered positive when flow is gravity-stabilized, whereas N b is considered negative when the gravity acting to destabilize the displacement. In the context of vertical gravity-driven DNAPL migration in porous media, N b is usually negative.

When two liquids are in contact in a porous medium, the denser fluid tends to sink to the bottom, and the upward force that this heavier fluid exerts on the lighter fluid to balance its weight is known as the buoyancy force. In the absence of buoyancy forces, the fluid patterns formed by immiscible displacement of a wetting fluid by a non-wetting fluid can be divided into three regimes (stable, viscous fingering, and capillary fingering regimes) separated by an intermediate transition zone. [START_REF] Lenormand | Mechanisms of displacement of one fluid by another fluid in a network of capillary ducts[END_REF] developed a phase diagram used to identify the location of the different regimes, and thus of the transitional region between them, in an M, N c space. [START_REF] Ewing | A generalized growth model for simulating initial migration of dense non-aqueous phase liquids[END_REF] noted that the phase diagram relating viscous, capillary, and stable flow regimes should contain a third dimension to account for gravity and buoyancy forces, given by the Bond number N b .

The quantitative parameters previously described have also been applied to estimate the mobility of DNAPLs in porous media in terms of temporal variation of DNAPL saturation profiles during immiscible displacements [START_REF] Pennell | Influence of Viscous and Buoyancy Forces on the Mobilization of Residual Tetrachloroethylene during Surfactant Flushing[END_REF][START_REF] Helen | Influence of viscous, gravitational, and capillary forces on DNAPL saturation[END_REF][START_REF] Hunt | Nonaqueous phase liquid transport and cleanup.I. Analysis of mechanisms[END_REF][START_REF] Theodoropoulou | Relative permeability and capillary pressure functions of porous media as related to the displacement growth pattern[END_REF][START_REF] Glen | Effects of initial saturation on properties modification and displacement of tetrachloroethene with aqueous isobutanol[END_REF]. Several methods have been used to study fluid saturations during the immiscible fluid displacement process. In the past, resistivity measurements [START_REF] Leverett | Flow of oil-water mixtures through unconsolidated sands[END_REF] and microwave absorption (Parson and Marathon, 1975) provided area-average fluid saturations. More recently, techniques such as radioisotope tracing/Gamma cameras [START_REF] Zhang | The velocity of DNAPL fingering in water-saturated porous media: laboratory experiments and a mobile-immobile-zone model[END_REF], X-ray tomography methods [START_REF] Glass | Gravity destabilized non-wetting phase invasion in macroheterogneous porous media: Experimental observations of invasion dynamics and scale analysis[END_REF] were able to deliver two-dimensional projections of the fluid saturation distribution and provide real-time visualization of the advancing front. In these cases, the ability to reproduce physical measurements was excellent. However, if excellent images were obtained and the exact macroscopic structure of fingers has to be provided, the application of this method was limited to either 2D flow cells, or fairly small 3D flow. One general conclusion from waterflooding investigations is that DNAPL saturation in porous media strongly depends on the interfacial tension between the immiscible liquids, the wettability of the solids with respect to the immiscible liquids, and the pore water velocity attained during water flooding [START_REF] Nguyen | The effect of displacement rate on imbibition relative permeability and residual saturation[END_REF].

The aforementioned methods have not provided enough accurate data for a conclusive study of displacement instability processes. In addition, few experimental studies on immiscible displacement have considered buoyancy forces, and relevant data from experiments depicting the development of DNAPL gravity fingering are thus needed.

Additionally, there is still no clear experimental evidence on the effect of flow rates, mode flow (gravity-stabilized, gravity-destabilized), and texture of the porous medium on the occurrence of this phenomenon. A consistent issue is that natural porous media are opaque, and we need to look at the pore scale to understand to understand how the fingers evolve.

The overall aim of the present work is to quantify displacement instabilities in a saturated homogeneous porous medium. The experimental approach is original as it is based on the use of optical fibers to quantify locally the arrival times of the DNAPL/water front at different points within a control section of the porous medium. The column experiments comprise a series of drainage experiments with TCE as the non-wetting fluid and water as the wetting fluid. A second set of experiments focus on a complete drainage-imbibition cycle to remobilize the TCE entrapped during the primary drainage in the porous medium. Various displacement configurations were considered for study of the influence of specific model parameters on the displacement process, such as the direction of flow, the velocity of the displacing fluid, and the texture of the porous medium. The experiments were analyzed with respect to both distribution of water/DNAPL front velocities, and the DNAPL saturation in the column. In addition, conclusions about displacement instability were then drawn.

Materials and methods

Experimental setup

Laboratory experiments were conducted on a transparent glass column with a total length of 68 cm and an inner diameter of 10 cm. Two porous media were used: medium sand and fine sand. The characteristic properties of both porous media, such as intrinsic permeability (k), mean grain size diameter (d 50 ), uniformity coefficient (U), and porosity (η) were determined experimentally (Table 1).

Trichloroethylene (TCE) was chosen as the DNAPL for the experiments because it is among the most frequently detected contaminants in subsurface environments. At 20°C, TCE has a density of 1.463 g/cm 3 , a viscosity of 0.0056 g/s/cm, and a solubility in water of 1300 mg/L [START_REF] Jellali | A large-scale experiment on mass transfer of trichloroethylene from the unsaturated zone of a sandy aquifer to its interfaces[END_REF]. In the column experiments, a dye (Organol Red BS) was added to TCE to allow visual detection of the DNAPL during the experiments; the dye concentration was approximately 5 mg/L. Organol Red BS is highly soluble in TCE, and does not modify its chemical and physical properties in the range of concentrations applied [START_REF] Jellali | Large scale experiment on transport of trichloroethylene in a controlled aquifer[END_REF]. A high-accuracy peristaltic pump with Teflon tubing and interchangeable heads for flow rates up to 200 mL/min was used to inject the displacing fluid in the sand-filled column. A 2-cm layer of coarse sand was placed at the top and bottom of the sand column to ensure uniform flow conditions in the inlet and outlet sections of the studied sand column. Pressures at the inlet and outlet were measured using a pressure transducer (Cerabar T -PMC131). 8 optical fiber sensors were developed (see Section 2.2) to measure the arrival time of the DNAPL front at different points within a section located halfway up the column. The measurement recordings were acquired using an Agilent 34970A Data Acquisition device connected to a personal computer. The experimental setup used for running the displacement experiments is illustrated in Figure 1.

The experimental setup was designed to produce a controlled and well-defined injection of TCE in an initially water-saturated homogeneous porous medium. Filling of the column was carried out with dry sand added regularly and successively in layers of several centimeters to a 10-cm deep water layer at the bottom of the sand column, which was then mixed from bottom to top with a thin rigid rod. Next, the sides of the glass column were gently tapped to encourage the sand to settle tightly. This method was previously demonstrated by [START_REF] Nsir | Etude expérimentale et numérique de la migration des polluants non miscibles à l'échelle de Darcy[END_REF] to yield packed columns with highly reproducible and uniform physical properties. To ensure that the sand filling had the same perfect water-wet characteristics in each run, the glass column was refilled with fresh sand for each experiment.

Before the start of each displacement experiment, the column was flushed with de-aired water to ensure complete water saturation of the sand, and to measure the hydraulic conductivity of the sand column. Once steady water flow was achieved, TCE was injected through the injection port of the packed column at a constant flow rate. The outlet section of the column was kept at a constant water pressure. The optical fiber sensors might disturb locally the flow field but the effect of their emplacement has no major impact on the conclusions drawn about displacement instability as their emplacements in the porous medium were not changed during the different studied displacement scenarios.

Two displacement modes were studied in the experiments, referred to as verticalupward and vertical-downward displacement later on in the study. In the first set of studies, drainage experiments were run in which TCE was the non-wetting fluid and water the wetting fluid. Here, the displacement is generally unstable according to the contrast of density and the unfavorable viscosity ratio M (M = 0.58<1) that existed between the immiscible fluids. The second set of experiments is primary drainage followed by water flooding to mobilize the TCE entrapped in the porous medium. The system constitutes a complete drainage-imbibition cycle. To reach residual TCE saturations, a large range of water flow rates was used. In practice, the critical value of the capillary number at which oil recovery becomes significant has been found to be in the range of between 10 -5 and 10 -4 . As the capillary number in our experiments is typically in the order of 10 -6 for ordinary water flooding, flow rate varied in the laboratory experiment from 40 up to 170 mL/min, corresponding to a capillary number N c of 4.1×10 -6 to 16.8×10 -6 . Residual TCE saturation after imbibition (water-flooding) was examined under a range of experimental conditions, including direction of flow (upward or downward), displacing fluid velocity, and the mean pore size (medium or fine sand). Table 2 and 3 summarize the different parameters varied in drainage experiments and experiments conducted to study drainage-imbibition cycle, respectively.

A major issue in the experimental studies was how to quantify the preferential pathways of DNAPL during its migration in a saturated porous medium. The fingering process was recorded using specially designed optical fiber sensors to measure the arrival time of the water/DNAPL front at eight points in a control section located halfway up the column. To obtain visual information on the fingering events that appeared in the sand column, digital images were taken of the red-colored DNAPL blobs observed at the glass wall of the column.

For each displacement experiment, the DNAPL breakthrough at the column outlet was determined by measuring the quantity of TCE leaving the column. The extent of local DNAPL saturation was determined by collecting and analyzing sand cores after each experiment (see Section 2.3).

Optical fiber sensors

The experiment used optical fibers to quantify the DNAPL migration pathways through the saturated sand. This technique, widely used in industrial applications, can act as an active sensing mechanism and is generally referred to as an intrinsic optical fiber sensor. Similar sensors have found application in chemical, biochemical, biomedical, and environmental sensing (Sharma and Guptal, 2006;[START_REF] Gupta | Surface Plasmon Resonance-Based Fiber Optic Sensors: Principle, Probe Designs, and Some Applications[END_REF]. The general structure of the developed optical fiber sensor is shown in Figure 2. It consists of an LED as the optical source, an optical fiber that transfers the optical signal, an optical detector (a photodiode), and processing electronics.

The sensing principle is based on the change of optical intensity modulation along the optical fiber. The tool converts the presence of a DNAPL into modified optical characteristics of the emitted signal. The original prototype version of the optical fiber sensor was conceived by VEGAS (Versuchseinrichtung für Grundwasser-und Altlastensanierung) of the University of Stuttgart [START_REF] Barczewski | Development and application of a light fiber fluorimeter for tracer tests[END_REF]. As part of the PhD thesis of [START_REF] Nsir | Etude expérimentale et numérique de la migration des polluants non miscibles à l'échelle de Darcy[END_REF], further improvements were undertaken to obtain a robust and sensitive tool for detecting the presence of DNAPL in porous media. The experimental work uses an all-silica multimode fiber with core/cladding/jacket external diameters of 600/675/690 μm. The optical fibers were approximately 50 cm long.

The sensing element was prepared by decladding a small portion of the optical fiber (about 1 cm in length), which was then cleaved at its end face. Next, the polished surface of the end face was coated with a thin aluminum layer to create a mirror, using monomer vapor phase deposition [START_REF] Nsir | Etude expérimentale et numérique de la migration des polluants non miscibles à l'échelle de Darcy[END_REF]. The in situ deposition of the chemically active metal (aluminum) on the polished surface of fiber is achieved by suspending it in the monomer vapor phase deposition. The output face of the fiber is placed in front of a crucible filled with boil-off and heated to above 1000°C. Evaporation of aluminum was carried out under a pressure of about 10 -7 atmospheres. The coating and adhesion of the metal to the fiber surface were of high quality.

The developed sensors do not indicate of what percent of the proportion of DNAPL was in contact with the sensitive part of the fiber. They react for primary droplets of DNAPL being in contact with the sensitive part of the fiber. Therefore, detection of DNAPL in any sensor emplacement does not mean that surface surrounding of the sensing element is completely impregnated by DNAPL. An estimate of DNAPL saturation is thus not possible with the actual fiber tool.

Quantification of DNAPL saturation

For each drainage experiment, a total volume of TCE corresponding to approximately 70% of the pore volume of the packed sand was injected into the column at a constant flow rate of 40 mL/min. The amount of TCE entrapped in the porous medium was determined by mass balance of DNAPL inflow and outflow.

Residual saturations of TCE were established in the water-saturated sand column after a complete drainage-imbibition cycle. First, TCE was injected into the saturated porous medium, as in the drainage experiment. The flow direction was then reversed, and water was pumped through the column to displace the mobile TCE phase. To minimize the removal of entrapped TCE by dissolution the water injected into the sand-filled column did not exceed three times the total pore volume. If no TCE mobilization was observed after flushing with the chosen flow rate, an increased flow rate was applied to the sand column. For each flow rate, the DNAPL recovered in the effluent was collected in glass vials, and the volume of TCE displaced as free product measured. The average TCE saturation in the column was calculated per difference.

At the end of each displacement experiment (drainage experiment/complete drainageimbibition cycle) the sand was sampled to obtain vertical profiles of local (entrapped/residual) TCE saturations in the saturated porous medium. Because of the emplacement of fibers, the sampling was achieved only above the sensor control section of the column. Local TCE saturations at various lateral positions and at different depths of the sand column were quantified using a soil sampler that probed a volume of roughly 2 cm 3 [START_REF] Nsir | Etude expérimentale et numérique de la migration des polluants non miscibles à l'échelle de Darcy[END_REF]. The sampling grid used for this process was very fine (5 cm spacing) in the vertical direction. At each sampling depth, one sample was taken on the axis of the column and two others on opposite points of the edges of the column section. TCE was extracted from these samples with methanol, and was measured by GC/FID.

Results and discussion

Arrival time of the water/DNAPL front for a low injection rate

To determine whether displacement instabilities occur, the arrival times of the DNAPL front were recorded at different points of the control section (z=-34 cm) during a drainage experiment conducted at a low injection rate of 40 mL/min, for both vertical-upward and vertical-downward flow modes. A normal distribution was fit to the arrival times, and normalized by the mean arrival time (Figure 3). All sensors placed in the control section of column detected DNAPL for both gravity-destabilized and gravity-stabilized flow. This may be explained by the fact that only a small portion of the sensitive part of fiber may have been in contact with DNAPL droplets. In the case of the downward flow mode, the arrival times of the DNAPL/water front varied widely, and their distribution is non-uniform. A standard deviation of nearly 70 sec was obtained, whereas in the case of upward flow mode a standard deviation of only 9 sec was calculated. This large difference is principally caused by densitydriven fingering that develops in the case of the vertical-downward displacement. The low injection rate makes viscous forces negligible, so the capillary and buoyancy forces together dictate the movement of the advancing fluid. Given the contrast of both density and viscosity between the two fluids, the injected non-wetting fluid (TCE) prefers to move along preferential paths with low capillary resistances, thereby forming pronounced fingers (Thibodeau et al., 1997). Capillary fingers, if formed, would be converted into gravity fingers as they elongated. A vertically oriented macroscopic finger may result that is composed of a sequence of larger cluster of entrapped TCE. In the case of vertical-upward displacement, gravity plays a stabilizing role in the migration of DNAPL and compensates for the destabilizing viscosity ratio M that exists between the two fluids. Under this condition, the water/TCE interface moves as a piston-like flow, and as a result, the recorded arrival times are very close together. This is consistent with the digital images taken of the red-colored TCE blobs observed at the glass wall of the column (figure 4). In the case of downward flow, TCE front was observed as small droplets or globules forming fingers and numerous protuberances, as described by [START_REF] Zhang | The velocity of DNAPL fingering in water-saturated porous media: laboratory experiments and a mobile-immobile-zone model[END_REF]. In the case of upward flow, the TCE front followed a compact path as a succession of homogeneous layers alternating with thin discontinuities. This layering is an artifact of the packing method employed in the experiments. Here, N b becomes positive in the system and its value of 2.7×10 -5 seems to be sufficient to generate a stable displacement regime independently of the values of both capillary number (N c =4.1 10 -6 ) and unfavorable viscosity ratio (M=0.058).

In case of downward displacement, buoyancy forces speed up the downward movement of the TCE on account of its higher density as compared to water (Figure 4b). Thus, the displacing fluid moves much faster than in the vertical-upward case, and a premature breakthrough time is established. This behavior was experimentally confirmed, with TCE breaking through at the outlet section of the column at nearly 30 minutes and 22 minutes for upward displacement and downward displacement, respectively.

Average DNAPL saturation

The average TCE saturations calculated by mass balance over the entire column are quite different for upward and downward displacements, and these observations confirm that density contrast plays a significant role in the migration of a DNAPL in a porous medium, particularly at low DNAPL injection rates. In the experiments, we determined an average TCE saturation of about 66% in the case of upward displacement, but only 33% for downward displacement. Furthermore, soil samples taken at the end of each drainage experiment clearly indicate that the TCE saturation front behavior is abrupt in the case of upward water displacement varying from 20% to 70% at a depth between 10 and 30 cm. In the case of downward displacement, the TCE saturations are more or less uniformly distributed over the sampling depths and are characterized by a low saturation of about 16% (see Figure 5). This difference is mainly due to the stable growth during upward displacement, whereas downward displacement is characterized by flow instabilities. The interactions between the capillary, buoyancy, and viscous forces result in a complex flow pattern characterized by ramified features as observed in figure 4b.

DNAPL pressure at the inlet section

Although the DNAPL saturation describes the distribution of displacement fluids in the porous medium, it is the pressure build-up that significantly determines the displacement of one fluid by another. Thus, to better understand the interplay between global pressure behavior and displacement condition (stable or unstable), the inlet TCE pressure was measured in the laboratory experiment for upward and downward displacement. Figures 6a and6b summarize the TCE pressure measured as a function of time at the inlet section of the soil column for the cases of upward displacement and downward displacement, respectively.

The boundary pressure applied to the outlet section of the experimental system corresponded to a water height of 18 cm and 86 cm, respectively. As shown in figure 8a, the pressure at the inlet section increased with time in the case of vertical upward displacement. Thus, the water/DNAPL interface is stabilized by the gravity effect in this case. This confirms the findings of [START_REF] Toussaint | Influence of porescale disorder on viscous fingering during drainage[END_REF]. The potentially destabilizing influence of the viscosity ratio was not significant. The pressure increase is therefore caused by the high density of the invading fluid, which overrides the decrease in viscous pressure drop during the displacement of the less viscous DNAPL. Downward displacement of water by TCE is an unstable displacement, with a much more irregular front and several gravity fingers. The interplay between hydrostatic pressure and viscous pressure varies under this condition of displacement and results in an overall decrease of the inlet TCE pressure. As shown in Figure 6b, the measured pressure decreases as the less viscous and denser fluid invades the system due to the constant TCE injection rate.

Here, buoyancy forces are destabilizing and dominate the movement, resulting in typical DNAPL fingers in the displaced fluid. The macroscopic water/DNAPL interface is therefore no longer horizontal and many throats become blocked by the capillary effects. The number of active invasion paths is therefore low, and the flow section of DNAPL is reduced. The contribution of viscous pressure is also limited and is dominated by the contribution of the higher density of the invading fluid. Combination of all of these effects results in the decreased inlet pressure.

Influence of the DNAPL injection rate on gravity fingering

The DNAPL injection rate is one of the key factors in displacement of water by DNAPL and controls the balance between the buoyancy, viscous, and capillary forces. Low injection rates yield more opportunities for capillary effects and buoyancy effects to acting in displacement, and the viscous forces are basically limited. It has been previously shown that under a low injection rate, the vertical-downward displacement of water by TCE gives rise to several fingers. Development of fingers is mainly caused by the destabilizing effect of the gravitational forces that dominate the displacement process in this situation. The invasion of a saturated porous medium by TCE is only partial, and the resulting distribution of the pollutant is thus inhomogeneous. To study the effect of increasing the injection rate on the occurrence of instability phenomena during downward flow mode, a displacement experiment was performed with a doubled DNAPL injection rate.

The statistical description of the measured arrival times shows that their distribution is quite similar to the one obtained for upward displacement at a low flow rate (see Figure 3).

The standard deviation is about 17 sec, which is significantly less than that obtained for the case of downward displacement of water by TCE at a low rate (standard deviation equal to 70 sec) and similar to that observed in the case of upward movement (stable case, standard deviation equal to 9 sec). Thus, the increase in the injection rate significantly stabilizes the displacement. The displacement pattern has more resemblance to the compounding pattern for the vertical-upward displacement recorded at low flow rate. This can be explained by the fact that the high DNAPL injection rate is related to a high-pressure gradient, and when the pressure gradient is high enough, and then both the capillary and buoyancy forces might be effectively neglected. Here N b is negative (-2.7×10 -5 ) and the unfavorable viscosity ratio M is higher, but an increase of the capillary number (N c increases from 4.1×10 -6 to 7.5×10 -6 ) seems to lead to a transition to a stable displacement regime which can be referenced in the 3D diagram proposed by [START_REF] Ewing | A generalized growth model for simulating initial migration of dense non-aqueous phase liquids[END_REF]. This result is rather surprising, because DNAPL has a lower viscosity, and one would expect a faster injection rate to destabilize the fronts even more. The motion may be unstable, but the shape of the finger head is such that it hits most of the sensors in a moderate amount of time, whereas the untouched parts of the flow are located in the zones without sensors.

Alternatively, the separation between a TCE finger and water may cover a longer distance at high speed, which then becomes comparable to or larger than the cell size.

Average TCE saturations obtained at high injection rate are shown as function of depth in Figure 5. The TCE saturations are globally higher than those obtained in the case of low injection rate. The corresponding profiles also display fewer irregularities than those obtained in the case of downward drainage conducted at a low injection rate. As the flow rate increases, the advancing DNAPL/water front behavior becomes more abrupt.

Effect of permeability on gravity fingering

To study the effect of the intrinsic permeability of the porous medium on the gravity fingering process, a second displacement experiment was performed using fine sand and applying the same low DNAPL injection rate as in the first experiment. The intrinsic permeability of the fine sand is approximately eighteen times lower than that of the medium sand. Figure 7 shows the comparison of measured dimensionless arrival times of the TCE front (during the downward displacement of water by TCE) carried out in fine sand to those obtained in the medium sand. In the case of the fine sand, statistical analysis of the determined distribution of front arrival times yields a nearly uniform distribution. The standard deviation in this case is only 12 sec, in comparison to the 70 sec recorded in the case of displacement of water by TCE in the medium sand. Because the permeability is primarily controlled by the pore-throat size, changing the size of the mean grain diameter of the sand has a larger effect on drainage displacement behavior due to the increase in capillarity effects [START_REF] Yao | Stability analysis of the unsaturated water flow equation 2. Experimental verification[END_REF]. One might expect that the growth of gravity-driven fingers would be nearly suppressed under this condition of displacement. Although buoyancy forces dominate the displacement especially at low DNAPL injection rates and the Bond number is negative, the capillary number (N c = 4.110 -6 ) is closer to the absolute value of Bond number (|N b |=4.210 -6 ) in the case of displacement in fine sand than in the case of displacement in medium sand ((|N b |=2.710 -5 ). Indeed, a lower mean grain size of the porous medium is directly translated to a lower absolute value of the Bond number (see Eq. 2). This might explain the reduction in the potential contribution of buoyancy forces to the development of fingering processes in low-permeability porous media. Furthermore, when an evolving finger propagates in fine sand, the finger velocity may be slowed down until the non-wetting phase pressure increases enough to overcome the local increase in capillarity effects, which is apparently sufficient to stabilize the displacement front. Indeed, the destabilizing effect of the unfavorable viscosity ratio on the displacement was not completely overcome and formation of local fingers that we did not observe might have occurred.

Based on these observations, we may assert that when the Bond number is negative but have the same magnitude as the capillary number, may corresponding to a stable displacement region that we can refine in the 3D diagram proposed by [START_REF] Ewing | A generalized growth model for simulating initial migration of dense non-aqueous phase liquids[END_REF].

Residual DNAPL saturations

Residual TCE saturations were measured in the sand-filled column after a complete drainage-imbibition cycle. Four configuration cycles related to the applied flow mode (vertical-upward flow, vertical-downward flow) in both primary drainage and water-flooding processes were carried out in the medium sand. The primary drainage was initially carried out under a stable displacement condition using a constant vertical-upward DNAPL flow rate. This resulted in a high initial TCE saturation, whereas the vertical-downward primary drainage gave rise to a small initial TCE saturation in the saturated porous medium (see Section 4.2). After each drainage experiment, water was then either injected from bottom to top, which reproduced an unstable gravity condition, or from top to bottom in a second experiment, allowing buoyancy forces to play a maximum stabilizing role in the water flooding process. In addition, one complete drainage-imbibition cycle was conducted in the fine sand (see table 3). In this case, both primary drainage and water-flooding displacements were performed in the direction of the gravity forces. Section 4.4 clearly showed that stable conditions were achieved during the drainage process even for a destabilized-gravity condition and unfavorable viscosity ratios between the two fluids.

A large range of water flow rates was applied to mobilize the DNAPL entrapped in the porous medium after primary drainage. Corresponding entrapped TCE saturation in the porous medium was measured by mass balance as described in Section 3.3. The water flow rate was varied (up to approximately four times the initial flow rate) by increasing the water injection rates from 40 to 170 mL/min, corresponding to a capillary number N c ranging from 4.1×10 -6 to 16.8 × 10 -6 . Figure 8 summarizes the observed average TCE saturations as a function of the capillary number (N c ) for the four drainage-imbibition cycles conducted in the medium sand and for the one cycle performed in both the fine and medium sand under the same flow mode.

It should be noted that, independent of the primary drainage condition (stable or unstable), the average residual TCE saturations quantified in the medium sand were lowest for downward water-flooding and highest for upward water-flooding at the lowest capillary number of 4.1×10 -6 . For a high initial TCE saturation in the porous medium (66%), the remaining TCE saturation is only 23 % for downward water-flooding against 45 % when a upward waterflooding was applied. This result is in good agreement with the findings of [START_REF] Helen | Influence of viscous, gravitational, and capillary forces on DNAPL saturation[END_REF]; In the case of upward water-flooding, the viscous forces are counteracted by buoyancy forces at the lowest water-flooding rates, and a large fraction of the TCE remains in the pores. However, buoyancy forces push the DNAPL droplets downward when the water is injected from the top. Whereas the buoyancy forces are cooperative forces in the downward water-flooding and are combined with the capillary and viscous forces to displace the DNAPL droplets, the buoyancy forces are subtractive forces in the case of upward displacement and compete the capillary and viscous forces that mobilize the entrapped DNAPL. Furthermore, the remaining TCE saturations in the medium sand for both flow modes applied during water-flooding tend to decrease as N c increases. This indicates that both the capillary and buoyancy forces are overcome by the viscous forces at high injection flow rates, and therefore, the mobilization of TCE becomes independent of the chosen flow mode.

On the other hand, it can be noted that residual TCE saturations obtained at high capillary numbers were slightly different for both water-flooding directions, whereas their values were fairly equal for the two primary drainage directions and followed the same water-flooding mode. Regardless of the initial saturation in the porous medium, the average residual TCE saturation was stabilized around 20% and 12% for upward water-flooding and downward water-flooding, respectively. Therefore, it can be concluded that residual TCE saturations are directly related to the imbibition flow mode and are independent of the initial distribution of DNAPL in the porous medium.

As shown in Figure 8b, the residual TCE saturation determined at the lowest capillary number (N c =4.1×10 -6 ) for the fine sand was slightly higher than that observed in the medium sand. The initial difference of about 5% increases with increasing capillary number. At higher capillary numbers, residual TCE saturation is reduced to 20% in the fine sand, and 13% in the medium sand. This may be explained by the fact that, with decreasing mean grain size of the porous medium, capillary forces become dominant and lead to a significant quantity of TCE being captured in the porous medium, even when viscous forces increase during waterflooding at high flow rates.

Figures 9a and9b show the measured local TCE saturations at a centre point and two point edges of a column transect as a function of depth for downward water-flooding displacement and upward water-flooding displacement in the medium sand, respectively. In the case of downward water-flooding, the residual TCE saturation seems to be more homogeneous along the whole section of the soil column because the measured TCE saturation values on the axis and edges (two opposite points) of the soil column are closer to each other at most sampling depths than those obtained in the case of upward-water flooding.

Despite favorable viscosity ratios existing between the immiscible fluids, fingering driven by destabilizing density may be formed [START_REF] Lenormand | Invasion percolation in an etched network: Measurement of a fractal dimension[END_REF][START_REF] Lenormand | Flow Through Porous media: Limits of Fractal Patterns[END_REF]. The shape of these fingers might be different to those observed during gravity destabilizing primary drainage. During downward water-flooding, as the displacing fluid (water) is lighter than the displaced fluid (TCE), the TCE/water interface is kept stable, and water most likely crosses the entire section of the sand-filled column. Consequently, the major part of the entrapped TCE in the porous medium is mobilized downwards at the same rate. However, due to the unstable gravity condition in the case of upward water-flooding, water flows upwards along preferential paths and leaves behind a heterogeneous spatial distribution of an amount of TCE in the sand. This is shown as an example for a depth of z = -15 cm (with z=0 corresponding to the top of the column). A low residual TCE saturation of only 2% was measured, indicating the preferential pathways of water, whereas at another sampling point located at the same depth, a high residual TCE saturation of nearly 25% was quantified, illustrating a region not in contact with the water flush.

The residual TCE saturation measured in the fine sand as a function of depth at different points of the column section (centre point and two point edges) are shown in Figure 9c. The residual TCE saturation appears to be more uniformly distributed than in the medium sand.

For a given depth, the difference between TCE saturation values in the axis and along the wall of the column is lower than 5%. Furthermore, the average residual TCE saturation varies between 16% and 24% along the total sampling height. This indicates that, due to the dominating capillary forces in the fine sand (as previously stated), the immiscible front of displacement is almost abrupt.

Conclusion

In this study, optical fiber sensors were developed to quantify DNAPL gravity-driven fingering in a water saturated homogeneous porous medium. The influence of several hydrodynamic parameters (such as direction of flow, velocity of displacing fluid, and texture of the porous medium) on the instability displacement process was evaluated during various displacement configurations (drainage experiment and complete drainage-imbibition cycle).

DNAPL front arrival times recorded halfway up in a column filled with medium sand were more spread out in downward drainage than those obtained in the case of upward drainage, showing evidence of preferential path growth for the non-wetting fluid in the saturated porous medium. Here, the spatial distribution of the measured DNAPL saturation was rather inhomogeneous. As expected, buoyancy forces help to stabilize the displacement during the vertical-upward drainage, however, they enhance the unstable displacement and may lead to gravity-driven fingers in the case of a downward displacement. It is worthwhile to note that the velocity of front displacement and the permeability of the porous medium may fundamentally affect the occurrence of gravity fingers and may even entirely suppress the fingering process during immiscible downward displacement. The experimental results highlight that, even for a gravity-destabilized condition, the displacement of the immiscible front was more stable in the case of a high DNAPL injection rate, and entrapped DNAPLs were nearly uniformly distributed in the saturated medium sand. This happens because viscous forces become dominant over buoyancy forces. This stabilizing effect is evident in laboratory experiments, but it is very hard to achieve in a real applied situation. Additionally, it was also concluded from the experimental results that downward migration of the DNAPL in a sand with fine texture leads to an enhanced stabilization of the immiscible front in the saturated porous medium. Here, the small pores of the fine sand accentuate the local capillary pressure in the pores, which slows down the front displacement and makes it more uniform.

The experimentally obtained residual DNAPL saturations clearly indicate that upward water-flooding displacement leads to higher quantities of entrapped DNAPL than during downward displacement of water by DNAPL. It was also found that residual DNAPL saturation is a decreasing function of the capillary number N c , which confirms the findings of earlier studies. Subsequently, the difference in permeability between the fine sand and the medium sand significantly impacts the values of residual DNAPL saturation. Furthermore, the residual DNAPL saturations measured in the case of the fine sand were much more uniformly distributed over the depth than in the medium sand. In an upcoming paper, the experimental data will be used to validate a recently developed numerical grain size distribution based pore-throat model [START_REF] Nsir | A pore-throat-model based on grain size distribution to quantify gravitydominated DNAPL instabilities in water saturated homogeneous porous medium[END_REF] to compute the water/DNAPL displacement at the pore scale.

Figure captions

  An optical fiber consists of a central core (through which light is guided) embedded in an outer cladding with a slightly lower refractive index. They are covered by a plastic jacket to mechanically protect the fiber and exclude stray light. Light transmission through an optical fiber occurs due to internal reflection at the core/cladding interface. Any change in reflection can be detected by a change in transmission. Because reflection depends on the refractive index of the medium surrounding the fiber, removing a portion of the cladding creates a zone that can detect DNAPL, due to the different refractive indices of water (n w =1.33) and DNAPL (n DNAPL = 1.46). Therefore, light intensity modulation can be brought about by refraction of rays in the DNAPL surrounding the modified cladding region.

Fig. 1 :

 1 Fig. 1: Schematic representation of the experimental setup.

Fig. 2 :

 2 Fig. 2: Principle of the developed optical fiber sensor. n cl and n c are the refractive indices of

Fig. 3 :

 3 Fig. 3: Measured dimensionless arrival times of DNAPL front in both vertical-upward flow

Fig. 4 :

 4 Fig. 4:TCE front observed on the column wall at the beginning (5 minutes) and the end of

Fig. 5 :

 5 Fig.5: Measured TCE saturation as function of depth and average TCE saturation in the

Fig. 6 :

 6 Fig. 6: Measured inlet and outlet pressure head as functions of time: (a) case of vertical-

Fig. 7 :

 7 Fig. 7: Measured dimensionless arrival times of the water/DNAPL front in both medium sand

Fig. 8 :

 8 Fig. 8: Average TCE saturations as a function of the capillary number N c : (a) four drainage-

Fig. 9 :

 9 Fig. 9: Residual TCE saturation at different positions of the column section as function of

Table 1 :

 1 Properties of the porous media used in the experiments

		intrinsic permeability	mean grain diameter	uniformity	porosity
		K [m 2 ]	d 50 [mm]	coefficient U [-]	(η)
	Medium sand	910 -11	0.43	2.1	0.43
	Fine sand	510 -12	0.17	2.3	0.40

Table 2 :

 2 Characteristics of drainage experiments

	experiment 1	porous medium medium sand	flow mode upward	N c 4.1×10 -6	M 0.58	N b 2.710 -5
	experiment 2	medium sand	downward	4.1×10 -6	0.58	-2.710 -5
	experiment 3	medium sand	downward	7.5×10 -6	0.58	-2.710 -5
	experiment 4	fine sand	downward	4.1×10 -6	0.58	-4.210 -6
	Table 3: Characteristics of drainage-imbibition cycle experiments
		porous medium	drainage	displacement	imbibition	displacement
			flow mode		condition	flow mode	condition
	experiment 1	medium sand	upward		stabilizing	downward	stabilizing
				gravity effect		gravity effect
	experiment 2	medium sand	upward		stabilizing	upward	destabilizing
				gravity effect		gravity effect
	experiment 3	medium sand	downward		destabilizing	downward	stabilizing
				gravity effect		gravity effect
	experiment 4	medium sand	downward		destabilizing	upward	destabilizing
				gravity effect		gravity effect
	experiment 5	fine sand	downward		destabilizing	downward	stabilizing
				gravity effect		gravity effect

Acknowledgements

Financial support for this research was received from the programme REALISE (REseau Alsace de Laboratoires en Ingénierie et Sciences pour l'Environnement). The Région Alsace, the GDR «Hydrodynamique et Transferts dans les Hydrosystèmes Souterrains » (INSU-CNRS), and the Conseil Scientifique de l'Université Louis Pasteur Strasbourg are gratefully acknowledged. Furthermore, the authors would like to thank the anonymous reviewers for their valuable comments and suggestions, which helped to improve the article significantly.