Alizée P Lehoux

Stéphane Rodts

Pamela Françoise Faure

Eric Michel

Denis Courtier-Murias

Philippe Coussot

Magnetic resonance imaging measurements evidence weak dispersion in homogeneous porous media

I. INTRODUCTION

In many environmental and industrial situations (ground water flows, chromatography, catalysis, etc) it is necessary to describe and predict the dispersion of solutes or contaminants that move with the flow of water through porous media. Mechanical dispersion, in which a solute or colloid particles in a fluid flowing through a porous medium progressively disperse while diffusion due to thermal agitation is negligible, is a well admitted concept [START_REF] Scheidegger | The physics of flow through porous media[END_REF][START_REF] Bear | Dynamics of fluids in porous media[END_REF]. Several theoretical approaches have been developed [START_REF] Bear | Dynamics of fluids in porous media[END_REF][START_REF] Hunt | [END_REF], mostly relying on the idea that since the structure is disordered each flowing element will follow a path similar to a random walk around the average velocity, a process leading on long times to Gaussian spreading. Some aspects of the physics of the process are nevertheless questioned [4] and anomalous or non-Fickian dispersion, generally attributed to porous medium heterogeneities, has been often observed [START_REF] Bear | Dynamics of fluids in porous media[END_REF]5] and is the subject of numerous theoretical developments [6]. It was even recently suggested that non-Fickian dispersion should be the norm [7].

Even for homogeneous systems, i.e. for which the physical properties of representative elementary volumes are homogeneous throughout the sample, the experimental knowledge appears somehow fragile. For similar conditions (material and flow characteristics) the basic parameter describing the process, namely the dispersion coefficient ( D ), is often not precisely determined [START_REF] Bear | Dynamics of fluids in porous media[END_REF]. The standard description and quantification of the phenomenon basically rely on the so-called "breakthrough curves", i.e. the time distribution, at the exit, of a solute injected as a step in the sample. This approach does not provide any information on the flow characteristics inside 3D porous medium. In this context existing data for flows through granular packings (of grain diameter d ) apparently fall, with some significant scattering, along a master curve [START_REF] Bear | Dynamics of fluids in porous media[END_REF]. In the regime of dominant mechanical dispersion this master curve can be represented by is the Peclet number, v the average (interstitial) fluid velocity through pores, and m D the particle diffusion coefficient, as estimated from particle size and Stokes-Einstein relation. Various original advanced approaches, relying on flow inversion experiments [8], particle tracking or scattering in matched-index medium [9], tracer concentration displacement [10] or imaging [11], provided further information on the solute distribution in time which apparently confirmed the above results.

In the other hand 1 H NMR (Nuclear Magnetic Resonance), which by means of PFG (Pulse Field Gradient) techniques can directly measure the statistics on molecular displacements in a flow over a given time interval, provides a powerful and precise technique to observe dispersion inside random bead packings. Generally, after travelling a distance of 10 to 20 bead diameters the displacement statistics in the main flow direction exhibits a Gaussian shape reminiscent of some asymptotic dispersion regime [12][13][14][15][16][START_REF] Stapf | NMR Imaging in Chemical Engineering[END_REF][START_REF] Akpa | [END_REF]. Apart from two tests [14], existing measurements with this technique generally provided values for smaller than that associated with the average master curve by a factor 2-3 [19][20], 4-5 [13], or 5-10 [21]. Such difference may be attributed to quite limited explored length scales, so that some macroscale effects responsible for higher macroscale dispersivity may be missed. Besides, from MRI (Magnetic Resonance Imaging) of D 2 O injection in chromatography gel columns a significant band broadening due to column inlet was observed [12], and from an analysis of NMR data, it was suggested that imperfect flow injection could also have a significant impact on the dispersion observed from breakthrough curve experiments [22]. This work also showed that a further analysis of PFG data makes it possible to distinguish intrinsic dispersion, which appears to be in agreement with the theoretical prediction for a random network of capillaries [23], whose resulting D is below the master curve (see above) by a factor about 3.

Thus there is a need to explore dispersion by direct measurements inside the sample over sufficiently large length scales and avoid potential impact of injection. In this work, we used MRI to get information on the internal flow characteristics and precise measurements of the evolution of the longitudinal concentration profile over time. We show that this profile is significantly affected by entrance effects which immediately deform the step radially. As a consequence a specific analysis taking into account this deformation must be set up to properly analyze dispersion. The new direct measurements finally obtained by this means provide effective dispersion coefficient values an order of magnitude lower than usually admitted.

II. MATERIALS AND METHODS

A. Materials

As porous media we used bead or sand packings. Glass beads supplied by Sigma Aldrich were washed twice with 5.8 M hydrochloric acid overnight and rinsed with pure water. Their mean diameter ( d ) was either 570 or 950 m. The size distribution (measured with a laser particle size analyzer) of such beads around these mean values are almost perfectly similar (see [24]) when rescaled by d , with a polydispersity coefficient, defined as

   i i i d d d P 
where i  is the volume fraction of the beads of diameter i d , about 10 % in each case. This means that these beads exhibit a slight polydispersity which precludes any crystalline arrangement in the packing.

The sand was "Fontainebleau sand" sieved between 200-250 µm, washed with heated nitric acid (65%) in a water-bath for 2 hours to dissolve organic matter and oxides, then with a 0.1 M NaOH solution to neutralize the acid, and finally rinsed with deionised water. Under these conditions the sand grain surface is negatively charged.

Sand or glass beads were carefully introduced in the PMMA column (length 21.5 cm, diameter 5 cm)) in order to avoid the presence of air bubbles. For this purpose, pure water was added from the bottom of the column in the same time as the porous material, with constant manual stirring of the mixture, in order to homogenize the sample and remove bubbles. Then the material was vibrated to enhance packing and finally the sample was closed with a cover pressed onto its top. This procedure leads to a packing for which relative displacements of grains are unlikely. The porosity of the sample in the column, computed from mass measurement of the different components (grains, water, empty column), ranged from 0.38 to 0.4 for the different samples.

As colloidal particles we used superparamagnetic commercially available particles (Molday ION Carboxyl) which have an iron oxide core and are functionalised with surface carboxyl groups (negatively charged). Their iron content is 1.49.10 -20 moles of iron per particle, and their density is 1.25 g/cm3. Diffusion Light Scattering (DLS) measurements at different concentrations showed that the particle diameter was constant and equal to 30 nm over 30 min., thus indicating that the particles don't aggregate. For these particles we have . The concentration of particles in water was measured from the iron concentration with D UV-visible spectroscopy Cary 50 Varian at 500 nm. These small superparamagnetic particles are homogeneously dispersed in water. They do not adsorb on the beads or on the sand grains. Moreover, due to a ratio between the grain size and the particle size between 7500 and 32000, most of particles are situated at very large distance (as compared to their size) from the pore walls and steric effects will certainly be negligible.

B. Experimental set up

The column containing the grain or bead packing is set up at the center of the magnet (see Figure 1). The column bottom is linked to a peristaltic pump. The liquid (pure water) or the particle suspension is injected from the bottom, and goes out of the sample at the top (see Figure 1), through a pipe (diameter 0 d =1.52 mm) along the central axis of the sample. The sample is ended at each end by a filter 420 µm thick with a pore size of 30 µm to avoid any loss of sand grain or glass bead through the pipe. During a steady flow of pure water through the porous medium packing we abruptly replace, for some time, water by a solution of colloidal particles at a concentration of 0.47 +/-0.04 Fe mmol/l at some distance in the upstream pipe (see Figure 1) then we resume the water injection. The injected volume of suspension, determined by weighting the beaker before and after the injection, varies slightly ( ml 2 . 0 

) around 10 ml depending on tests. This volume is equal to 6% of the total liquid volume in the column. The characteristics of the tests are presented in Table 1.

Due to the flow along the pipe length (2 m) between the point of injection and the sample entrance there might be some effect of the upstream flow on the (initial) distribution of particle concentration reaching in the sample. Since the flow is laminar (   0

Re

Vd  is equal to 70 for a typical velocity of 5 cm/s, where V the mean velocity in the pipe) the flow is essentially a Poiseuille flow with a parabolic profile (

  ) / 1 ( 13 2 R r 
in m). This implies that the fluid arriving at the column will be a two-phase flow with a central part made of suspension and an outer part made of pure liquid. During the time needed to reach the entrance there will be a limited exchange between these two phases. Indeed the typical length of diffusion over a time of 40 s is 17 microns, which is almost one hundred times smaller than the pipe diameter. After the entrance in the column the flow characteristics are complex and we could expect some kind of mixing of these two phases leading to an apparent concentration smaller than the injected concentration. It is likely that there is a mixture of these two phases inside the porous medium at the scale of pore. However, if this effect had an impact on the concentration profile it should be more important for a longer distance of pipe between the injection point and the entrance, but comparing the two series of concentration profiles (inside the sample) obtained for a distance (pipe length between injection point and entrance) of 2 and 5 m we do not see any difference. A possible explanation is that this effect remains minor as compared to the strong entrance effects (flow enlargement, see main text). Under these conditions we assume that, in the same way, the flow in the pipe after the exit has no impact on the results in terms of breakthrough curve.

For some of these tests we also measured the so-called breakthrough curve. A fraction collector was set up at some distance (2.9 m) from the exit (see Figure 1) which allows getting fluid samples. Output were sampled regularly (approximately every 8 ml) and weighted. The iron concentration of each sample was measured from light extinction in visible spectrophotometry at 500 nm, and the breakthrough curves could be drawn from these data.

C. MRI

MRI experiments were carried out with a vertical imaging spectrometer DBX 24/80 Bruker operating at 0.5T (20MHz proton frequency) and equipped with a birdcage radio frequency coil of 20 cm inner diameter. The suspended colloidal particles modify local NMR relaxation times of water molecules depending on their local concentration, and can then be detected or quantified by means of relaxation weighted 1 H MRI [25]. 2D Images were taken using 2D NMR slice-selective spin-echo imaging at 6.54ms

 E T , with 300ms  R T
. The chosen field of view was 24*10 cm, the matrix size 256*64 pixels, and the slice thickness 10 mm. Fastness constraint of the 2D imaging procedure prevents variations of signal intensity in such pictures from being strictly proportional to particle concentration. Quantitative measurements can nevertheless be achieved on 1D "profiles" MRI pictures, where the sample is projected onto its longitudinal axis x . The local intensity



. Note that instead of a simple spin echo, this MRI protocol uses the second echo of a CPMG (Carr-Purcell-Meiboom-Gill) sequence [26] to better compensate field inhomogeneities. In the above equation, 0 S is the signal intensity of the slice without particles, R the particle relaxivity calibrated from independent tests with particle suspensions, and ms 41 .

3

 E T
is the CPMG echo time. Here the matrix size of 128 pixels and the recycling delay

s 0 1  R T . Short E T allows approximating   E RT x S x S x 2 / ) ( / ) ( ln ) ( 0   
. Even in the limit of short E T , the latter formula was indeed found far more reliable than a just a linear approximation in the case of inhomogeneous tracer repartition through the slice. Prior to data analysis, axis distortions in x direction at bottom and top of the sample are calibrated and corrected. Measurements of   ) (x  over long experimental times also suffer baseline fluctuations attributed to temperature evolution in NMR coils and electronics on the time scale of few hours. Note that these fluctuations were taken into account and corrected in our analysis (see Section III.B). To sum up, from such measurements, we get the profile of concentration distribution along the sample axis, i.e. the average particle concentration in successive crosssectional layers along the axis x .

III. EXPERIMENTAL RESULTS AND DISCUSSION

A. 2D Imaging

2D MRI in a central longitudinal section of thickness 10 mm provides successive pictures of the spatial particle distribution along the sample axis (see Figure 2). Around the entrance and the exit the images are deformed due to magnetic field heterogeneity but a relevant qualitative analysis of the flow characteristics is still possible: the first images show that the front of particles is not straight. We could check that this shape does not result from the upstream flow in the pipe as it was insensitive to the distance between the sample entrance and the injection point in the pipe (see Section II.B).

Actually this front shape finds a simple explanation: as the fluid initially arrives from the small pipe and enters the porous medium along its central axis ( x ) the fluid elements which will invade the sample sides have to move radially and thus will initially advance in the longitudinal direction more slowly than fluid elements around the central axis. This effect may even yield dead zones in the corners of the entrance [27], which does not seem to occur in our case. Such effect is inherent to injection in a porous medium and a previous attempt [10] to sort it out by setting the solute as a straight step at the sample bottom was probably unsuccessful because of deformation occurring during flow start. It also seems that any type of grid or chamber positioned just before the entrance in the porous medium cannot solve this problem. Indeed, as long as the flow remains laminar and particle diffusion is negligible at the pore scale, in steady state, the streamlines, whatever their shape, will need to cover larger distances to reach the lateral edges, so that in general a similar effect as that observed here should occur. The only way to rule out this problem would consist in pushing the liquid, and then the suspension uniformly distributed in a cross-sectional layer, with a piston of the same size as that of the sample. However this seems very difficult to achieve in practice. At last note that the front shape as shown in Figure 2 is the general case (see also Figure 3) but it can also be more complex as observed with sand in the two series of MRI pictures taken with this material (see Figure 4). In that case we suspect that the sample could exhibit some slight packing heterogeneity around the entrance which allowed an apparently better initial radial spreading but we have no more precise explanation. Anyway these observations do not affect the analysis and results described below. Finally, after a displacement over a distance approximately equal to the sample diameter the radial shape of the particle distribution in the 2D layer apparently negligibly evolves (see Figures 2,3,4): it essentially spreads longitudinally but its curvature remains constant. A more precise analysis confirms this statement. The front shape, simply taken as the apparent outer limit of the yellow regions, was determined from the successive pictures such as those shown in Figures 2,3,4. The corresponding curves are plotted in Figure 5 for data of Figure 2. The curves superimpose very well within the noise on measurement, it is not possible to detect any trend of evolution of the global shape of the profile, whatever the flow direction and the position of the pulse in the sample. These observations mean that on average: i) the particles essentially move longitudinally;

ii) their average velocity in this direction is now uniform;

iii) the longitudinal dispersion, which was already at work during the entrance stage, is now the major effect.

These conclusions are confirmed by further tests. If we now reverse the flow direction just before the particle front reaches the sample top we do not see any further significant evolution of the radius of curvature of the front (see Figure 2,3,4 d,e,f,g) and this is the same for the next flow reversal (see Figures 2,3,4 g,h,i,j). Thus, as long as the step does not get out of the sample there is no further front deformation but only dispersion. Moreover it is remarkable that there is apparently no impact of gravity and no impact of the flow reversal (indeed the dispersion coefficient was independent of the flow direction, see below). This shows that the dispersion in that regime and over a length much larger than the grain size is fully irreversible, confirming the previous observation of Rigord et al. [28] from ionic tracer echo dispersion measurement. Since the transverse dispersion coefficient was observed, either by standard techniques or by NMR [13], to be more than ten times smaller than the longitudinal dispersion coefficient, we will neglect it here. Actually this assumption is consistent with our 2D observations, in particular the fact that the front shape remains apparently constant when moving three times through the sample (see Figures 2,3,4 then Figure 5) whereas significant transverse dispersion with such a curved pulse front (associated with a transverse particle concentration gradient) would tend to progressively deform the front.

There is nevertheless an effect which deserves to be mentioned. In the pictures (especially Figures 2 and3) we can see a spike along the wall, essentially visible on the "left" side, whose height grows as the particle pulse advances. It has been suggested that the dispersion is enhanced due to radial heterogeneities in bead packing such as those resulting from the wall proximity, even for rather large ratios between the column diameter and the bead diameter [29][30]. However in our case this ratio is of the order of 1/100, so that this effect should be minor. Moreover here this spike has two remarkable characteristics: (i) its growth rate is larger than that of the pulse width, so that it reaches a larger height than the central pulse front at the approach of the column end (see Figures 2d and3d); (ii) this spike tends to disappear during the reversal flow (see Figure 2,3 e,f,g), and then again significantly grows during the last upward flow (see Figure 2,3 h,i,j). This means that, in contrast with dispersion (see above), this effect is almost fully reversible and tends to develop more widely in the flow direction. Such an effect is likely due to some inhomogeneity of the sample along the wall, inducing a slightly smaller solid volume fraction and a larger permeability allowing locally higher velocity. Nevertheless, due to the limited resolution of these 2D pictures, it is not relevant to try to further quantify this effect from our data. As a consequence we will not take it into account in the data treatment below but we can keep in mind that this effect tends to slightly enhance the apparent dispersion, and thus give slightly larger values for the dispersion coefficient than for a flow through a perfectly homogeneous sample.

As already mentioned, fastness constraint of the 2D imaging procedure prevents the quantitative conversion of signal intensity in such pictures to particle concentration. We can nevertheless get a qualitative view of the concentration distribution along the sample axis by integrating the signal, obtained in these 2D images, in each elementary cross-section. The resulting total signal associated with each position gives a profile of the concentration distribution along the longitudinal direction. From the successive pictures along the main axis we get a series of profiles (see Figure 6) whose height slightly decreases while their width increases. Dispersion is thus already apparent.

B. 1D imaging

Concentration profiles

From the 1D imaging procedure we can get a quantitative description of the evolution of the particle distribution along the sample axis. Typical results (see Figure 6) show a good consistency between these profiles and those deduced from imaging (which as for them were extracted from 2D images). Differences owe to shifted timings for 1D and 2D measurements, lack of quantitativity of 2D data, and the fact that unlike 1D profiles, 2D measurements do not probe the whole sample section. This confirms that the shape of the 1D profiles is essentially determined by the step deformation at the entrance and then dispersion can tend to further spread them. of paramagnetic volume particle concentration (  ) at different times during the first upward flow (every 89 s) for the test of Figure 1. Profiles (dashed red lines) computed from 2D images (see text) (every 106 s). For these profiles the concentration scale (assumed to be proportional to the grey level) was arbitrarily chosen so as to get a level for the initial profile close to that from 1D profile at the same instant.

Qualitatively similar results are obtained for different materials (sand packing and glass bead packing) with mean grain diameters from 220 to 950 µm, different step volumes and average fluid velocities from . This corresponds to Peclet ( Pe ) numbers between 200 and 20000 (see Table 1), a range in which mechanical dispersion is dominant [START_REF] Bear | Dynamics of fluids in porous media[END_REF]. In order to illustrate the range of aspects of concentration profile evolutions during the flow through the column due to the fluctuations from one test to another, we show three additional typical examples, corresponding to approximately the same Peclet number but different resulting values of dispersion coefficient (see Table 1). To better show the profile shape evolution we shifted them of the displacement expected from the mean flow velocity recorded (see Figure 7). x (cm)

) mol.l (Fe -1 c -4 -2 0 2 0.0 1.0x10 -4 2.0x10 -4 (c) x (cm)

Analysis of dispersion

As the 2D images do not provide full quantitative information on the dispersion inside the whole sample we will rely on the 1D profiles. However, since these profiles are generally strongly impacted by the heterogeneity of the 2D distribution of colloids, the analysis of dispersion from such profiles has to take this into account. In that aim we assume that initially the particle concentration in water is of no particular shape due to uncontrolled entrance effects. Regions occupied by the suspension are supposed to be distributed in the form of longitudinal fingers of some constant characteristic thickness and progressing at the same velocity v along sample axis. In the sample let f s be the cross section area of a given finger f , and ) (x F f the local tracer concentration in this finger at 0  t (see Figure 8). x direction, the new profile at time t writes:

  ) , ( ) ( ) , ( ) ( ) , ( 0 t x g x F s t x g x F t x F f f f      (1)
where  is a convolution product along x direction, and

Dt Dt vt x t x g  4 / ) 4 / ) ( exp( ) , ( 2    is the 1D dispersion propagator. ) , ( t x F
then still obeys the 1D equations for Gaussian dispersion. The dispersion coefficient D is determined from a least square fit to the whole data set, where the initial concentration profile ) ( 0 x F , v and D are taken as fitting parameters. Note that the important parameters sought for in this analysis are transport parameters D and v only, ) ( 0 x F is deduced from the fitting procedure and does not need to be measured directly, and a precise definition of time origin is not required, provided it is taken after the pulse enters the sample, and before the first profile is measured. Note that we could check that the fitted velocity was close to the value measured in the hydraulic network. The fitting procedure is described in detail in Appendix 1.

This approach applied to each of our tests well fitted the data within experimental noise level, and provides the longitudinal dispersion coefficient with a maximum uncertainty of 10% (error bars covered by symbols in Figure 9). Crossed symbols correspond to tests with a sequence of upstream-downstream-upstream flows with the same particle step. The grey area represents the region covered by previous data (by conventional techniques) as gathered in Figure 1 of [16], which covers the data leading to the "universal" flow curve as for example shown in [START_REF] Scheidegger | The physics of flow through porous media[END_REF] and various other data obtained more recently. Dispersion coefficients determined from breakthrough curves are represented by filled symbols (single flow) and half-filled symbols (two flow reversals). The continuous line is equation [START_REF] Bear | Dynamics of fluids in porous media[END_REF].

The values for the dispersion coefficient obtained by this approach for all our tests are shown in Figure 9. No specific difference is observed for the different material types (grain size, narrow or wide grain size distribution): in logarithmic scale the dispersion coefficient falls along a straight line with a significant scattering. We have no clear explanation for this scattering: it is usual to consider that it is due to fluctuations in the detailed characteristics of the sample structure from one sample to another; however, depending on experiments, sequences of tests with successive upstream and downstream flows through the same sample may exhibit either no difference at all in the dispersion coefficient (superimposed crossed symbols) or some significant variation from one flow to another of the order of that which may be observed for different samples (see Figure 9), suggesting that other effects occur. These effects might be due to slight velocity heterogeneities in a crosssection not taken into account in our modeling.

Under these conditions it is very instructive to estimate the dispersion coefficient via the standard procedure, i.e. from the breakthrough curve. We find values which are significantly scattered and, on average, situated slightly below those generally shown in literature (see Figure 9). These values obtained from breakthrough curves are finally significantly larger than those deduced from internal measurements by a factor of about 4. This result is consistent with our observation that the shape of the concentration profile is mainly governed by entrance effects. Thus the breakthrough curve results from this initial deformation, then some dispersion in the sample, and finally a possible further deformation at the exit since the fluid elements the farthest from the central axis will reach the exit later than those situated close to this axis. As a consequence the impact of entrance and exit effects can a priori be minimized by using a sufficiently long sample. This is indeed what we observe from measurements of dispersion coefficient from breakthrough curves after two flow reversals, and thus a longer path, inside the column: the resulting coefficient values are smaller and thus closer to the values determined from internal observations (see Fig. 9).

Finally the dispersion coefficients obtained by a direct analysis of the evolution of the concentration distribution inside the sample, and a priori not affected by any artefact, are significantly smaller (on average by a factor about 10) than values obtained so far from standard experiments (see Fig. 9). As a first approximation we can represent these data by the following simple function:

Pe D D m 4 . 0  (2)
This tends to confirm the lower than previously admitted dispersion coefficient deduced from NMR PFG experiments at the length scale of a few tens of glass beads (see Section I). Note that it has been suggested that the Schmidt number (i.e. ratio of kinematic viscosity to diffusion coefficient) has an impact on dispersion [31][32]. According to these references the dispersion coefficient is larger for higher Schmidt number, but the correlation suggested in this publication does not go beyond a Schmidt number of 2000. If we nevertheless extrapolate qualitatively this result to our situation which corresponds to a Schmidt number of 54000, we can conclude that our experiments should provide larger dispersion coefficients than usually found, which tend to reinforce our conclusion on the discrepancy between our data and usual data obtained from breakthrough curves.

IV. Conclusion

By directly monitoring the motion of nano-particles over the sample length, and analyzing for the first time these displacements with a straightforward method insensitive to entrance or exit effects, we find dispersion coefficient on average lower than the data usually serving as reference for engineering (e.g. models predicting pollutant dispersion in environment) by almost an order of magnitude. Previous NMR-PFG measurements over relatively short distances inside the sample also provided lower values of the dispersion coefficient, but in a range slightly above our average value. A possibility is that many previous data based on the standard interpretation of breakthrough curve were affected by entrance and exit effects. This questions the relevance of previous theories attempting to represent these former data. Moreover the fact that dispersion directly measured at much lower values than usually expected in homogeneous systems suggests that it might be a rather difficult task to appreciate the effective "abnormality" of dispersion in other systems.

Acknowledgments:

We are grateful to J.P. Hulin for enlightening discussions, and we thank Y. Méheust, T. Le Borgne and R. Juanes for fruitful discussions during the initiation of this work. We acknowledge technical help from B. Jouaud, G. Sevenier and M. Vielpeau from EMMAH. where Shape(x) can be modelled as either a 6th order polynomial, or a Fourier expansion truncated to its 7 first terms.

Appendix 1. Fitting formulas for 1D profile analysis

The sum of were used to ensure fitting to be well posed and stable. The fitting algorithm was an iterative constrained trust-region method, and the iteration was stopped when least square criterion could not be further reduced regarding computer double precision accuracy. It was always checked visually that fits well corresponded to data. Residual fitting discrepancy was found to well correspond to noise level on measurements, and this noise level was then used to estimate error bars on fitted parameters. In particular, relative errors on dispersion coefficient never exceeded 10%. It was finally checked that fitted D values did not significantly depend on base-line description, nor on the detailed choice of

Figure 1 :

 1 Figure 1: Scheme of principle of the experiment.

Figure 2 :

 2 Figure 2: (color online) Images of the particle distribution along longitudinal cross-sections at successive times during the flow through a bead packing (570 µm): first upward flow (a, b, c, d) ( -1 4 m.s 10 14 . 1    v

Figure 3 :

 3 Figure 3: (color online) Images of the particle distribution along longitudinal cross-sections at successive times during the flow through a bead packing (950 µm): upward flow (a,b,c,d), downward flow (e,f,g), final upward flow (h,I,j). Same scale and colour scale as in Figure 2. First upward flow (a, b, c, d) ( -1 4 m.s 10 06 . 1    v

Figure 4 :

 4 Figure 4: (color online) Images of the particle distribution along longitudinal cross-sections at successive times during the flow through a sand packing (225 µm): upward flow (a,b,c,d), downward flow (e,f,g), final upward flow (h,i,j). Same colour scale as in Figure 2. First upward flow (a, b, c, d) ( -1 4 m.s 10 05 . 1    v ), downward flow (e, f, g) ( -1 4 m.s 10 21 . 1    v

Figure 5 :

 5 Figure 5: (color online) Successive shapes of the front during pulse displacement in the sample for the test shown in Figure 2: first upwards motion (red continuous), downward motion (blue dashed), second upward motion (black dotted).

Figure 6 :

 6 Figure 6: (color online) Longitudinal 1D profiles (continuous blue lines, corrected from baseline fluctuations)

Figure 7 :

 7 Figure 7: Shifted (see text) concentration profiles for different tests with close Peclet number values: (a) large beads with

Figure 8 :.

 8 Figure 8: Scheme of the principle of calculation described in the text.

Figure 9 :

 9 Figure9: Dispersion coefficient scaled by the diffusion coefficient of the particles computed from the 1D NMR profiles for the different tests as a function of the Peclet number for the different materials: sand (squares), small beads (circles), large beads (diamonds). Crossed symbols correspond to tests with a sequence of upstream-downstream-upstream flows with the same particle step. The grey area represents the region covered by previous data (by conventional techniques) as gathered in Figure1of[16], which covers the data leading to the "universal" flow curve as for example shown in[START_REF] Scheidegger | The physics of flow through porous media[END_REF] and various other data obtained more recently. Dispersion coefficients determined from breakthrough curves are represented by filled symbols (single flow) and half-filled symbols (two flow reversals). The continuous line is equation[START_REF] Bear | Dynamics of fluids in porous media[END_REF].

  a sum of neighbouring positive Gaussian functions following: list of about 15 fixed coordinates spanning the initial peak domain with a constant step dx . We used cm 0.5  dx in all our fits. It was checked that smaller values did not improve fitting procedure, and that changing the interval covered by n x 's did not produce any change in deduced transport coefficients. With such a Gaussian decomposition, model function for ) , ( t x F is still analytic, following: to their specific thermal nature, baseline distortions in our measurements for each measuring time i t are modelled as the product of a fixed baseline shape and a fluctuating coefficient:

nx

  's -provided they well span the initial peak domain-, nor on small variations of constraints put on n a 's and n  's, i.e. variations were kept far under the error bar.

3122 Table 1 :

 31221 Test parameters and measurements. The relative uncertainty on measured dispersion coefficient is +/-10% (see text and appendix).

  Shape. The fit usually involved about 2000 to 40000 data points, and 100 fitting coefficients. Relevance of the fit was then carefully controlled. First of all, constraints

				F	, ( t x	)	and	, ( t x Base	)	are fitted to the data, taking as fitting parameters	D , v	,	a 's,	 's ,	p 's,
													n	n	i
	and polynomial or Fourier coefficients of	(x	)
													a		0	and
													n
	3 dx	 	n 	dx	/	3


T Chevalier

C Chevalier

X Clain

J C Dupla

J Canou

S Rodts

P Coussot

Darcy's law for yield stress fluid flowing through a porous medium

We measured the pressure drop vs flow rate during the flow, in a wide range of velocities, of well controlled yield stress fluids through confined packings of glass beads of different sizes. A detailed analysis of the data makes it possible to extract a general expression for the pressure drop vs flow rate curve through a porous medium as a function of the flow rate and the characteristics of the system. This general law has a form similar to the Herschel-Bulkley model describing the rheological behavior of such fluids in simple shear, i.e. it expresses as the sum of a critical (yielding) pressure drop and a flow rate dependent term. This law involves the rheological parameters of the fluid, one characteristic length of the medium, and two coefficients which only depend on the structure of the porous medium. The first coefficient is related to the path of maximum width throughout the porous medium while the second coefficient reflects the pore size distribution. The values of these coefficients were determined in the case of a granular packing.

Introduction

The flow of non-Newtonian fluids through porous media is of interest for various applications: penetration of glue in the surface porosity of solid materials, injection of muds, slurries or cement grouts to reinforce soils, propagation of blood through kidney and, likely the most economically important application, injection of drilling fluids in rocks either for the reinforcement of the wells or for enhancing oil recovery, etc. Much work has focused on the flow properties of purely viscous non-Newtonian fluids, the simplest example being power-law fluids [START_REF] Chhabra | Flow of non-Newtonian fluids in fixed and fluidised beds[END_REF][START_REF] Liu | On non-Newtonian fluid flow in ducts and porous media[END_REF][START_REF] Lopez | Predictive network modelling of single-phase non-Newtonian flow in porous media[END_REF][START_REF] Sabiri | Pressure drop in non-Newtonian purely viscous fluid flow through porous media[END_REF]. For example the propagation of a non-Newtonian liquid pushing another one, and the possible instability of the interface leading to fingering, is a subject of importance which has already been studied [START_REF] Chen | The flow and displacement in porous media of fluids with yield stress[END_REF][START_REF] Zitha | Control of flow through porous media using polymer gels[END_REF]. Another aspect concerns the effect of inertia which can further complicate the analysis [START_REF] Morais | Non-Newtonian fluid flow through three-dimensional disordered porous media[END_REF]. However there is a wide range of injected fluids, and in particular in oil industry, which exhibit a strongly non-Newtonian behavior, i.e. they are yield stress fluids which flow like liquid only when the stress applied to them is larger than a critical value. These fluids are emulsions, foams, clayey suspensions or even some polymeric gels. As far as the flow of such fluids through porous media is concerned this yield stress character is likely to induce specific properties which are so far poorly identified. The basic knowledge of the pressure drop needed to induce a steady inertia-less flow of a yield stress fluid through a porous medium is so far not very well developed. In this paper we mainly focus on this situation. In practice the ultimate goal will be the possibility to predict the pressure drop vs flow rate relation from independent measurements of the rheological characteristics of the fluid and of the geometrical characteristics of the porous media (based on macroscopic parameters such as a permeability).

The existing approaches for establishing the pressure drop vs flow rate relation for yield stress fluids have relied on the usual physical approaches of the permeability coefficient appearing in the Darcy's law for Newtonian fluids: they use the analogy of the flow at a local scale within the porous medium with the flow either through a capillary or around isolated particles. Thus one can obtain an equation relating the pressure drop to the average velocity involving specific but partly unknown lengthscales (pore size, total flow path length). In this way various expressions for ) (V f p   have been obtained [START_REF] Al-Fariss | Flow through porous media of a shear-thinning liquid with yield stress[END_REF][START_REF] Chase | A correlation for yield stress fluid flow through packed beds[END_REF][START_REF] Park | The flow of non-Newtonian fluids through porous media[END_REF][START_REF] Pascal | Nonsteady flow of non-Newtonian fluids through a porous medium[END_REF][START_REF] Vradis | Macroscopic conductivities for flow of Bingham plastics in porous media[END_REF] which all have the following structure, typical of yield stress fluids: the velocity differs from zero only beyond a critical pressure drop c p  ; then it progressively increases as the difference between the actual pressure drop and c p  increases.

The relevance of these approaches for yield stress fluids is questionable. In particular Balhoff and Thompson [START_REF] Balhoff | Modeling the steady flow of yield stress fluids in packed beds[END_REF] suggested that for each of these relationships "it is not yet clear whether universally applicable adjustements (in particular of the lengthscales) can be made or whether a case-by-case assessment is needed". Balhoff and Thompson [START_REF] Balhoff | Modeling the steady flow of yield stress fluids in packed beds[END_REF] first evoke that the converging-diverging behavior of individual throats has a prominent effect for yield stress fluids.

In fact we emphasize that the essential point is that in contrast with the flow through a capillary, the flow of a yield stress fluid through a pore does not necessarily involve a constant size of the possible volume of fluid transported: as the pressure drop increases a wider region of fluid starts to flow. So, in contrast with Newtonian fluids, as the pressure drop increases the characteristic length related to the fluid volume available for flow at a local scale would increase. A second difference of flow characteristics is that at a macroscopic scale flow actually starts as a percolation effect: at the critical pressure drop the flow occurs only along a specific path throughout the porous medium [START_REF] Balhoff | Modeling the steady flow of yield stress fluids in packed beds[END_REF][START_REF] Sahimi | Non-linear transport processes in disordered media[END_REF]; then, we can expect that as the pressure drop is increased more flowing paths progressively form within the porous medium. In fact these effects should not lead to an apparent behavior at a macroscopic scale fundamentally different from that which has been observed or predicted: they are consistent with a critical pressure threshold and then an increase of pressure drop with velocity. The point is that the lengthscales involved in this law might not vary at all as expected in the above theoretical approaches.

Besides there have been several approaches based on network models [START_REF] Chen | The flow and displacement in porous media of fluids with yield stress[END_REF][START_REF] Sochi | Pore-scale network modelling of Ellis and Herschel-Bulkley fluids[END_REF] which provide more straightforward descriptions of the reality but don't provide generic analytical expressions for the pressure drop vs velocity relationship.

Unfortunately experimental studies in that field [START_REF] Al-Fariss | Flow through porous media of a shear-thinning liquid with yield stress[END_REF][START_REF] Chase | A correlation for yield stress fluid flow through packed beds[END_REF][START_REF] Park | The flow of non-Newtonian fluids through porous media[END_REF] are still scarce and rely on parameter variations (pore size or rheological parameters) in rather limited ranges. In addition existing data provide measurements of flow through porous medium in rather narrow ranges of mean velocity, typically slightly more than one decade. This implies that the comparison with a model can hardly be discriminatory. Also the comparison is done in linear scale so that the discrepancies at low flow rates are neglected, a critical point which precludes a relevant estimation of the critical pressure drop. At last the authors don't discuss whether the ranges of shear rates explored in both flow types (rheometry and flow through porous medium) are similar so as to support the relevancy of a comparison with the predictions of a model using the rheological parameters measured independently.

Actually a major difficulty in experiments of non-Newtonian fluid flow through porous media is to be able to obtain reliable data bases for both the rheological behavior of the fluid and its transport properties through the porous medium. As it is well known rheometrical tests with yield stress fluids need to be carried out carefully in order to obtain a relevant description of the material behavior [START_REF] Coussot | Rheometry of pastes, suspensions and granular materials[END_REF] but on the basis of current knowledge it is possible to get reasonably relevant data. With an experiment of yield stress fluid flow transport through a porous medium various additional problems have to be taken into account. For example: standard pressure measurements cannot necessarily be used because of the development of normal stresses in cavities; additional pressure drop in the conduits before and after the porous medium may be significant for a yield stress fluid; the elasticity of the fluid in its solid regime may induce transient effects; the rheological properties of the fluid may be affected by its contact with the solid surfaces of the porous medium (ions or impurities) (a classical effect with Carbopol gels).

We propose here an attempt to gather a set of relevant data for the two types of tests (rheometry and flow transport through porous media), for a wide range of shear rates and different sizes of porous media, different pore sizes and fluid yield stresses. With that aim, we used a set up in which the different possible experimental problems mentioned above are minimized. The wide data base obtained under such different conditions makes it possible to analyze the results in detail and get a general expression for the pressure drop vs flow rate, i.e. an empirical Darcy's law for yield stress fluids, with the appropriate scalings. After presenting the experimental set up, materials used and procedures in Section 2, we then present the data in Section 3, which are subsequently analyzed in Section 4.

Materials and procedures

Materials

In order to validate our measurement system, procedure and data analysis approach we carried out tests with a (Newtonian) Glycerol solution (99+%, Alfa Aesar), exhibiting a viscosity of 1.5 Pa.s at 20°C.

The first type of model yield stress fluid used was Carbopol (U980) solution in water. It has been observed that this material is essentially a glass made of a high concentration of individual, elastic sponges (with a typical element size of 2 to 20m) [START_REF] Piau | Carbopol gels: elastoviscoplastic and slippery glasses made of individual swollen sponges meso-and macroscopic properties, constitutive equations and scaling laws[END_REF], which gives rise to its yielding behavior.

The key advantage of this material is that it shows limited thixotropy [START_REF] Piau | Carbopol gels: elastoviscoplastic and slippery glasses made of individual swollen sponges meso-and macroscopic properties, constitutive equations and scaling laws[END_REF][START_REF] Coussot | Aging and solid or liquid behavior in pastes[END_REF] and phase migration, making it a quite well-behaved yield stress fluid. The preparation of Carbopol gel begins with the introduction of some water in a mortar mixer. The rotation velocity is set at 90 rpm and the appropriate amount of raw Carbopol powder is slowly added to the stirring water. After about one hour, the incorporation of the powder is done and the appropriate amount of Sodium Hydroxide (1 mol/L) is quickly added to the solution, which increases its pH. The mixing is then maintained for approximately one day to allow a full homogenization of the mixture. Such a procedure makes it possible to obtain a material which remains stable for several months. Different Carbopol concentrations were used, leading to different rheological properties. Concentrations from 0.1 to 0.8% were used but because of a significant influence of the actual mixing procedure, mixer type and mixing speed, on the solution finally obtained, we do not use, in the following, the concentrations as a parameter to identify a given solution but we rather rely on its rheological properties.

The second type of model yield stress fluid was water-in-oil emulsion obtained by dispersing in an oil phase droplets of brine stabilized by a surfactant. At high volume fraction of the dispersed phase, droplets come into contact and the fluid exhibits a yield stress [START_REF] Mason | Yielding and Flow of Monodisperse Emulsions[END_REF]. Batches of water-in-oil emulsion are prepared by dispersing a water + CaCl 2 (150 g/l) salt solution in a surfactant solution (dodecane (Acros Organics) + 7.5%wt of Span 80 surfactant) with a Silverson L4RT mixer. The rotation speed of the mixer is kept at 700 rpm during the addition of brine. It is then increased up to 6000 rpm during 85 min until a homogeneous emulsion is obtained. The droplet concentration was varied slightly around 80%. In that case a slight evolution of the system (possibly some kind of coarsening) could be observed between some periods of tests. As a consequence, like for the Carbopol solutions these materials were basically characterized through their rheological properties. The droplet size was around 1m.

Rheology

Rheological tests were performed with a Bohlin-stress-controlled rheometer equipped with two circular, rough plates (diameter: 40 mm). The sample was carefully set up and the gap was fixed at 2 mm taking care not to entrap air bubbles. A logarithmically increasing then decreasing stress ramp test was then applied over a total time of four minutes. Except for the first part of the increasing curve associated with deformations in the solid regime the increasing and decreasing shear stress vs shear rate curves almost superimpose. We retained, here, the decreasing part as the flow curve of the material. For similar Carbopol solutions and emulsions it has been shown that this apparent flow curve obtained from macroscopic observations correspond to the effective, local constitutive equation observed at a local scale with imaging technique [START_REF] Coussot | Macroscopic vs. local rheology of yield stress fluids[END_REF].

Carbopol gels and emulsions exhibit a simple yield stress fluid behaviour and their flow curve can be well fitted by a HB (Herschel-Bulkley) model (see typical results in Figure 1):

n c c k          (1)
in which  is the shear stress, 0   the shear rate, c  the yield stress, k the consistency factor and n the power law exponent. For the different Carbopol solutions the material parameters were found in the following ranges:

Experimental set up

The setup used for carrying out the experiments has been adapted from a set of experiments initially developed for studying the injection of micro-cement suspensions into sands, taking into account the filtration and clogging phenomena [START_REF] Dupla | An advanced experimental set-up for studying monodirectional grout injection process[END_REF][START_REF] Saada | Evaluation of elementary filtration properties of a cement grout injected in a sand[END_REF][START_REF] Saada | Modelling of cement suspension flow in granular porous media[END_REF], with application to the in situ improvement of mechanical properties of soils using grouting/injection processes. The set up is made up of two main components, (see Figure 2) which are: (i) an injection setup composed of a two-part (oil/fluid) injection reservoir, connected to an air/oil interface cell allowing to apply a regulated oil pressure in the reservoir; (ii) a confinement cell which contains the confined porous medium. This cell is connected, through plastic tubings, to the injection reservoir. The experimental setup is instrumented with pressure transducers (measurement of oil an fluid pressures during the test) and with a precision weigh scale installed at the exit of the confinement cell (see figure 2.a) for continuously recording the mass of fluid flowing out of the porous medium and calculate the flow rate.

The injection reservoir (see Figure 2a) is composed of two parts (bottom and top parts) separated by a low stiffness waterproof "Bellofram" membrane. The reservoir is 15 cm in diameter and 36 cm in height. The bottom part of the reservoir is filled up with oil (Newtonian fluid). This part is connected to the air-oil interface cell, which makes it possible to control the injection pressure. A pressure transducer is placed at the oil outlet point. The latter measures the relative oil pressure, p i , with respect to the atmospheric pressure. The increase of oil volume in the bottom part of the reservoir induces a displacement of the membrane which, in turn, pushes the fluid to be injected, located in the upper part of the reservoir, out of the reservoir. The fluid flowing out is then channeled to the confinement cell and injected in the porous medium. We believe that the present technique used for measuring the injection pressure, via the upstream pressure in a Newtonian fluid (the oil), makes it possible to avoid possible errors which could occur if placing a pressure transducer directly in contact with the yield stress fluid. A simplified cross-section of the confinement cell, showing the porous sample, is presented in Figure 2b. The porous medium is contained in a thin waterproof latex membrane which allows the application to the specimen, through the water filling up the confinement cell, a confining pressure Pc. Metallic grids are placed at the top and bottom of the specimen in order to retain the granular material (beads) during the test. The diameter of the specimen is 8 cm, different heights may be used (4, 8 or 16 cm). In order to prevent the development of preferential paths and obtain a uniform flow entering the specimen, the latter is placed on a conical hollow bottom plate. A similar conical end plate is adjusted on top of the specimen for the outgoing fluid.

The fluid flowing out of the sample falls into a small reservoir set on the precision balance mentioned before (Figure 2.a). We checked that the tested fluid did not evolve during the test by sampling it at the exit tube and comparing its rheological properties with those of the injected fluid.

Porous medium

We used different porous media made up of glass bead packings. Different types of beads were used, with a mean diameter ( D ): ranging from 0.26 to 2 mm. We measured the size distribution of the beads with a laser granulometer (Laser Beckman-Coulter). Typical sets of data obtained with this technique are shown in Figure 3. It appears that the size distributions of the beads are almost perfectly similar. A more detailed analysis of these distributions can be carried out with the help of the polydispersity coefficient defined as

    i i i D d D P   100
where i  is the volume fraction of the beads having a diameter i d . We found that P varies between 8 and 12% depending on the grain type. This confirms that our grains are almost monodisperse with a very similar polydispersity, which justifies that in the following we consider that we have a single type of porous medium, i.e. the different porous media are similar and characterized by only one length scale, namely the mean diameter D . As some impurities remaining on the bead surface may destabilize Carbopol solutions we tested the impact of cleaning the beads in an ultrasonic bath during one hour. The data obtained with or without such cleaning were similar.

The porous medium is made by slightly compacting, manually, successive layers of beads in a forming mould against which is initially applied the latex membrane. The porosity finally obtained is % 6 34 . 0  , which in particular means that our porous medium preparation was reproducible. Then, a confining pressure kPa 100

 c P
is applied to the porous medium, which moulds the membrane onto the beads and maintains the grain packing. It is expected that this technique limits the development of preferential flow paths at the periphery of the porous medium, a phenomenon which tends to occur when the bead pack is in contact with a planar surface. Indeed, in the latter case, the porosity is larger close to the planar surface than inside the porous medium. Note that we did not observe any significant compaction or widening of the sample section due to the confinement.

In order to calibrate the system for head losses due to the equipement itself (calibration for system compliance) and to be able to make corresponding corrections on the raw data obtained, specific tests were carried out by replacing the porous medium by a rigid empty tube that was called the phantom. The phantom, equipped, inside, with a latex membrane similar to the one used for the real test, had the same inner diameter as the porous medium. The cylinder was then externally covered with another latex membrane, allowing to adjust the phantom with the bottom and top end plates of the confinement cell.

Experimental procedure

We first slowly introduce the fluid in the top part of the injection reservoir. For yield stress fluids this step is crucial in order to limit the number of air bubbles introduced into the reservoir during this operation. Then the porous medium is prepared as described above. The injection reservoir is then connected to the confinement cell. Afterwards the impregnation of the porous medium can start: the oil pressure is increased to reach a pressure gradient of approximately 20 times a critical pressure gradient (estimated from the fluid yield stress divided by the bead's mean diameter). Such a pressure a priori makes it possible to saturate the porous medium in a short time with a sufficiently slow flow rate. The process is stopped when the flow observed on the exit tube of the confinement cell is well established. At that time the sample is supposed to be fully saturated. Note that at any time during the test it is necessary to check that the confining pressure is significantly larger than the oil pressure in order to maintain the rigidity of the bead packing.

We then start to increase the pressure step by step until reaching a continuous flow at the exit. Then the pressure is decreased step by step. This decrease is stopped when the out-flow does not seem to stabilize, typically for a flow rate less than 0.3 drop/min. For each pressure level, the measurement is initiated after stabilisation of the flow as observed on the exit tube, and this pressure is maintained during a sufficient time to get a measurable mass evolution at the exit. For example, for flow rates as low as a few drops per minute we wait for at least 10 drops going out of the device. Note that a different bead pack was used for a test with each new fluid, i.e. after a series of measurements at different velocities with one fluid type the beads were thrown away.

Experimental results

Procedure for analyzing the data

A typical data set for p i (t) and m(t) is shown in Figure 4. These data can be analyzed in terms of pressure drop per unit length, i.e.

h p p i  
in which h is the height of the porous medium, vs the flow rate, i.e. the measured (average) velocity of the material through the porous medium:

t m S V     1 ( 2
)
in which S is the section area of the porous medium and  the fluid density. Note that the above definition for the pressure drop will find a full physical meaning only when we have computed the effective pressure drop associated with the flow through the porous sample only (see below). A typical apparent flow curve, i.e. p  vs V , obtained in this way is shown in Figure 5. ) including a pressure (circles) increase step by step followed by a pressure decrease, and the corresponding mass evolution in time (squares). For each injection test with beads we also carry out a calibration test with the phantom, and analyzed the data the same way. We thus obtain two apparent flow curves (see Figure 5). In each case the pressure drop includes different pressure drops due to the flow through the different parts of the system except the porous medium. It also includes the hydrostatic pressure drop. These contributions are a priori the same when one uses a porous medium or a phantom. These calibration corrections allow us to obtain the effective flow curve, i.e. strictly associated with the flow through the porous medium (and without hydrostatic contribution), by subtracting the apparent flow curve associated with the phantom data from the flow curve found with the porous sample. Such a procedure in particular makes it possible to get relevant data without having to measure the pressure in the yield stress fluid (for example at the inlet and outlet of the confinement cell). This operation would indeed be difficult and not necessarily reliable because the pressure is not well defined in a yield stress fluid at rest. Note that in practice in order to facilitate the subtraction we first fitted a model to the flow curve data for the phantom (see Figure 5). . Estimating the typical shear rate in the porous medium from the ratio of the effective mean velocity through the porosity to the typical pore size (around one third the bead radius) we find that in all our tests it varied in the range This corresponds to the range of shear rates covered in our rheometrical tests.

We observed flow rates in the range

Results for Glycerol

Injection tests were carried out with Glycerol for two bead sizes (1 mm and 2mm). After subtracting the (constant) hydrostatic pressure term the flow curve obtained for the phantom appears to be linear (see Figure 6). This is the typical result obtained for the flow through a conduit or a porous medium. This indicates that in our range of measurements the whole system around the bead packing behaves like a given porous medium, which supports our procedure for extracting the effective flow curve (see above). The effective flow curves obtained for the Glycerol solutions through the porous media are also linear (see Figure 6), in agreement with Darcy's law:

V K P    (3)
in which K is the permeability of the medium. Since the porous media are similar (see Section 2.4) except for their characteristic length ( D ) we deduce from scaling arguments that the permeability should be proportional to 2 D :

2 D K   .
From our data (linear fits of Figure 6) we find for  . This globally confirms the validity of our approach but shows that the results are not perfect, which appears reasonable considering the uncertainties at the different steps of the process (phantom, slight viscosity variations, slight differences in the porous media, etc). Note that the empirical Kozeny-Carman equation [START_REF] Kozeny | Über Kapillare Leitung des Wassers im Boden, Sitzungsber[END_REF] for the same porosity predicts a slightly smaller value ( -4 x10 5 



) than the above one, which may be due to some difference in the grain size distribution.

It is worth recalling that for a liquid flowing through a straight cylindrical conduit the flow curve is formally similar to (3) but now with



. This means that the equivalent diameter of a conduit which would represent the porous medium in terms of Darcy's law for a Newtonian fluid is equal to 6

. 6 ' D D d    
. The value of the final factor, i.e. 6.6, between the bead diameter of the porous medium and the diameter of the equivalent conduit in terms of Darcy's law mainly results from the three following effects: (i) an average thickness of the paths between the beads significantly smaller than the bead diameter; (ii) a section area available for fluid smaller than the apparent section due to the presence of the beads, which increases the effective velocity through the equivalent conduit; (iii) the tortuosity of the fluid paths which tends to increase the length of these paths.

Results for yield stress fluids

The apparent flow curves for the phantom have a shape similar to the flow curve of the material obtained from rheometrical tests. Under these conditions we decided to fit these curves with a Herschel-Bulkley model with the same exponent as that used for fitting the flow curve (i.e.

36 . 0  n ).
The resulting effective flow curves for the different materials and mean bead diameters all present the same shape: they exhibit a plateau at low velocities and the pressure drop progressively increases with increasing velocity. This shape is similar to the typical (rheometrical) flow curve of a simple yield stress fluid. This means that in order to induce a flow through a porous medium one needs to overcome a critical pressure drop. When the bead size decreases or when the material yield stress increases the shape of the flow curve does not seem to change significantly (see Figures 7 and8): at first sight it seems that this is essentially the critical pressure drop which increases while the additional velocity dependent component is approximately constant (the curves are roughly similar by vertical shift in a logarithmic scale). The more detailed analysis presented below will show that this is not always right: the velocity-dependent component also significantly varies when the pore size varies. Anyway the vertical shift (in logarithmic scale) of the flow curves for a decreasing pore size or an increasing yield stress globally corresponds to the expected impact of the permeability decrease or the viscosity increase as for any type of viscous fluid flowing through any type of porous medium. Let us now compare the results obtained for Carbopol solutions and emulsion having similar rheological properties and flowing through the same porous medium. Surprisingly the effective flow curves are significantly different (see Figure 9). In particular the critical pressure drop for the emulsion is significantly lower than for that for the Carbopol solution by a factor of about two. This difference cannot be explained solely by the difference of values of the parameter k for the two types of fluid. This suggests that there is an unexpected impact of the nature of the fluid, a point that we will discuss further below. Considering the shape of the effective flow curves for all the injection tests it is quite natural to fit a Herschel-Bulkley model to the data keeping the same exponent as in the fit on rheometrical data. This model here expresses as:

n p c V k P P     (4)
in which c P  is the critical pressure drop and p k the consistency parameter for injection. This step is important as it means that it will now be possible to fully analyze the results simply through the values of the two parameters of this model as a function of the experimental parameters (bead diameter, fluid yield stress and consistency, flow rate).

In a first step we can now plot the two injection parameters as a function of the corresponding rheological fluid parameters. Globally each parameter increases when the other one increases but it is clear that they is no straightforward exclusive relationship between them. The data representing the pressure drop as a function of the yield stress for all our experimental conditions are widely dispersed.

Discussion

Dimensional analysis

A simple comparison of the HB model parameters for injection and rheometrical tests is not relevant. Indeed c P  and p k also depend on the characteristics of the porous medium through which the fluid flows. In order to take this aspect into account we can again take advantage of the fact that the porous medium is characterized by a single lengthscale, i.e. D . From a dimensional analysis we readily deduce that we should expect two relations in the form:

D P c c     (5) 1   n p D k k  (6)
where  and  are two unknown parameters which are independent of the material type and porous medium characteristics.

In order to check the validity of this analysis we plotted the data in the form suggested by equations ( 5) and [START_REF] Zitha | Control of flow through porous media using polymer gels[END_REF]. From this representation we can conclude that

c P  is indeed proportional to D c  and p k is proportional to 1  n D k
(see Figure 10): all the data obtained for each material fall along a straight line of slope 1 on a logarithmic scale. For p k a single straight line is obtained for the data corresponding to both materials, but for c P  two different straight lines are obtained. This is consistent with our above observation, i.e. the injection flow curves obtained for an emulsion and for a Carbopol solution with similar yield stress were significantly different. Finally our results can be summarized in terms of a general model giving the pressure drop as a function of the flow rate, i.e. the equivalent of the Darcy's law for yield stress fluids:

n c D V k p D            (7)
In this expression  and  should be universal factors for a porous medium composed of an assembly of spheres, and D V is the apparent shear rate for the flow through such a porous medium. The first (yielding) term on the right hand side of equation ( 7) corresponds to the critical pressure for incipient flow. As usual the second term, which is velocity dependent, derives from the additional viscous stress above the yield stress as long as the flow has been initiated.

The differences observed for the values of  for the two types of fluids would imply that they do not see the same porous medium, which is obviously not the case. Therefore, from a general point of view we need to understand the origin of the values for  and  and it is expected that this understanding will provide some elements of explanation for the discrepancy found between the two fluids. 4)) fitted to data of injection tests as a function of the appropriate expressions deduced from dimensional analysis and involving the corresponding rheological parameters of the fluids for Carbopol solutions (squares) and emulsions (stars) and the porous medium characteristics. The continuous lines correspond to straight lines of slope  equal to 12 for the Carbopol and 5.5 for the emulsion in (a), and of slope 85   in (b).

Analogy with the flow through a simple conduit

In order to further analyze our results it is useful to compare them with the flow characteristics for a flow through a straight conduit. Let us consider the steady uniform flow of a Herschel-Bulkley fluid in a cylinder of radius R . The momentum equation can easily be solved to obtain the mean velocity (flow discharge per unit surface) as the function of the pressure gradient (see for example [START_REF] Bird | The rheology and flow of viscoplastic materials[END_REF]). This may be expressed in terms of the relationship between two dimensionless numbers, namely

c p R G  2  
and

n n c kV R Bi   :                          ) 3 )( 2 )( 1 ( 2 ) 1 ( 2 ) 2 )( 1 ( ) 1 ( 1 ; 0 1 2 1 3 1 m m m G m G m m G G Bi G Bi G m m (8) in which n m 1  and L p p   
where L is the conduit length. An explicit analytical expression for the pressure gradient as a function of the flow rate cannot be found from this equation. However we can deduce the expressions in the asymptotic cases:

c R p Bi  2      (9a) n n R V k R m p Bi            1 ) 3 ( 2 0 (9b)
Note that the expression (9a) could be obtained directly from a simple momentum balance for describing the incipient motion of the yield stress fluid through the conduit. On another hand the expression (9b) is simply the solution of the flow problem with a power-law fluid. ) : critical pressure for incipient motion (dotted line), viscous term in the absence of yielding (dashed line), sum of the two previous terms (thin continuous line). The exact theory (according to equation ( 8)) is represented by the thick continuous line.

An approximate expression can naturally be found by adding the two asymptotic terms above and we obtain:

n c d V k p d          ' '    (10) in which R d 2  , 4 '  and   n n 2 6 4 '   
. This expression appears to well represent the theoretical curve within 10 % over the whole range of Bi (see Figure 11).

It is instructive to compare the values obtained for Darcy's law for the porous medium to those corresponding to the flow through this simple conduit. With [START_REF] Park | The flow of non-Newtonian fluids through porous media[END_REF] we have an expression equivalent to that obtained from our experiments for the flow through a porous medium so that we can find directly the conduit size which would lead to the same pressure drop vs flow rate expression. Let us first consider the yielding term, i.e. c  . For the Carbopol solution and the emulsion we find values (see Figure 10a) which lead to a ratio between the bead and the conduit diameters which are respectively 1.4 and 3. These values significantly differ from the factor (6.6) corresponding to the Newtonian Darcy's law (see Section 3.2), which suggest that there is no simple analogy between the physical mechanisms in the two situations.

Let us now further examine the possibility of equivalence between the two types of flow. This may be done by remarking that for the flow through a simple conduit the ratio of the two factors is constant, i.e. independent of velocity, rheological parameters and conduit size: for emulsions. This clearly means that, although it gives reasonable values for the critical pressure drop, the equivalence between the simple conduit flow and a porous medium for a yield stress fluid is not consistent.

Analogy with the flow through a bundle of conduits

In order to understand the origin of the above observed inconsistency we look at the flow characteristics of a yield stress fluid through a bundle of parallel conduits of different sizes. We consider a set of

N conduit of widths R i  ) ,... 1 ( N i 
. The number of conduits per size ( i n ) is inversely proportional to their section area so that the total area for each conduit type is constant. All the conduits are parallel to each other, so that the applied pressure drop is given and independent of the conduit size. Under such conditions, as we increase the pressure drop the flow is initiated in the conduits of largest size while the fluid remains in its solid regime elsewhere. Then as the pressure drop is increased the flow begins in the slightly smaller conduit and so on.

We define a specific dimensionless pressure drop for each conduit radius:

G i G i   ) (
, in which G is the value for the largest conduit. The flow in each conduit type is governed by the value of

) (i G according to: 1 1 ) ( ; 0 1 ) (       i i x i G x i G
. From the mean velocity through each conduit (8) we deduce the mean velocity through the set of conduits:

                                 ) 3 )( 2 )( 1 ( 2 ) 1 ( 2 ) 2 )( 1 ( ) 1 ( 1 ; 0 1 2 2 1 2 3 1 m m m G m G m m G x n n G Bi G Bi G i i i m i i i i i i m     (11)

in which

Bi is defined the same way as above, based on the radius of the largest conduit. Finally the pressure drop vs velocity curve has the usual form with a plateau associated with the critical pressure drop for flow to occur in the first conduit and a finite slope for large velocities. This slope should tend to 1 for very large velocities when the fluid flows in all of the conduits. As a consequence we can still use the approximation of the sum of the two asymptotic expressions for the pressure drop (as in [START_REF] Park | The flow of non-Newtonian fluids through porous media[END_REF]):

m i i m i i i i m n n Bi m G / 1 3 2 1 / 1 ) 3 ( 1                (12)
With such a form the characteristics of the set of conduits is accounted for by a single coefficient multiplied by

1 
Bi . This approximation is slightly not as good as for the single conduit, the maximum error being 15% for 36 . 0  n in the transition region, but still constitutes a useful relevant representation. Finally this gives an expression equivalent to [START_REF] Park | The flow of non-Newtonian fluids through porous media[END_REF] for the pressure drop with now:

4 '  and m i i m i i i i m n n m / 1 3 2 / 1 ) 3 ( 4 '               (13)
The origin of this new value for '  is that an increasing amount of fluid flows when the pressure drop increases. We computed the value of the ratio ' '   for a range of one hundred conduit radii distributed in a logarithmic way between the maximum value d and the minimum value d  (see Figure 12). We see that widening the range of diameters induces, indeed, a decrease of the ratio ' '   . This is not sufficient to reach a value as low as that observed for the emulsion (0.065) but it appears possible to reach an appropriate value by using another arbitrary type of size distribution, in particular by taking into account a higher fraction of lower diameter conduits.

Finally we can conclude that the first coefficient in the pressure drop equation, i.e. the equivalent of Darcy's law for yield stress fluid, mainly finds it physical origin in the widest path through the medium, while the second coefficient finds its origin in the complex multiscale flow through the medium. As a consequence the basic origin of the relatively low value of the coefficient  is the fact that in a porous medium the widest fluid path is first mobilized for the critical pressure drop while thinner and thinner paths, and thus smaller and smaller pores are then mobilized for increasing pressure drop. The next step will be to discuss the origin of the different values found for  . The curious point is that we have a consistency of data through the general scaling laws with regards to a continuum approach, whereas these laws differ for two materials with different internal structures but the same mechanical properties. 12) as a function of the minimum to maximum diameter ratio for the flow through a series of conduits of different diameters (see text).

Flow through sieves

In order to understand the origin of the discrepancy between the values of  associated with the critical pressure drop for the two different materials we carried out tests consisting to leave the material flow through a simpler porous medium, i.e. a sieve, now under the action of gravity only. The experimental device is composed of a plexiglass tube (30 cm long and 8 cm wide) over a sieve with a mesh size ranging from 20 m  to 2 mm. The set-up is weighed during the flow in order to obtain the mass of fluid percolating through the sieve.

Figure 13 shows a comparison of the results obtained for a Carbopol solution with a yield stress of 65 Pa, and an emulsion with a yield stress of 76 Pa, for two different sieve sizes. The mass evolves by jumps because the fluid going through the sieve progressively accumulates just below it and ultimately falls when it has reached a critical mass value. This resembles the behavior of a simple liquid falling by droplets but here the origin of the process is a competition between gravity and yield stress, which leads to larger "droplets". In order to clarify the representation of data we only plotted one symbol per level reached after each new droplet drop except when the duration between two successive drops was too long (in that case we plotted several points, which generally corresponds to an apparent stoppage of the flow).

We can see that despite a slightly higher yield stress, the flow of Carbopol through a given sieve is faster and stops before the flow of the emulsion through the same sieve. This means that the Carbopol solution needs a higher pressure drop to flow than the emulsion. Actually this effect is almost negligible for large mesh sizes (say, beyond m 250  (see Figure 13)) and increases when the mesh size decreases, becoming dramatic for a mesh size below m 20  . In the latter case there is almost no flow observed through the sieve for the Carbopol while the flow is still significant for the emulsion.

These experiments show that even for a simpler geometry there is a difference between the flow characteristics observed for the two fluids as soon as the pore size becomes sufficiently small. It is very likely that this effect be at the origin of the differences observed between both fluids since, roughly speaking, the critical mesh size can be considered similar to that of the smallest pores. Thus there might be some kind of jamming of Carbopol blobs in the smallest pores, an effect which would occur for smaller mesh sizes for the emulsion made of smaller elements.

However the point which remains very surprising is that the ratio between the values for  for the two fluids apparently remains constant over a wide range of pore size (related to bead diameter). If this was solely an effect of jamming it would be significant for small pores and would become negligible for sufficiently large pores. This means that further studies are needed to completely understand the origin of this effect.

Conclusion

We have carried out flow tests through well-calibrated porous media (glass beads), in which both the rheological behavior of the fluids and the flow characteristics through the porous medium (beads packing) were well-controlled. A thorough analysis of the data made it possible to extract the general pressure drop vs flow rate curve for a yield stress fluid flowing through a porous medium. This general law contains a yielding term which may be simply expressed as a function of the yield stress of the material and the bead size. The results obtained with the emulsion seem more reliable than those obtained with Carbopol solutions, which were likely to be affected by some jamming effects in the smallest pores. As for a Herschel-Bulkley model, an additional viscous term, depending on the flow rate, had to be added to the yielding term in order to obtain the total pressure drop during flow, which gives the general Darcy's law for yield stress fluid flowing through a bead packing:

n c D V k p D          85 5 .

 (14)

This expression is expected to be valid only if the elements composing the yield stress fluid are much smaller than the pore size, which obviously can never be perfectly true.

The yielding term in ( 14) is related to the mobilization of the fluid in the widest path throughout the porous medium. As the flow rate increases it is expected that additional thinner fluid paths progressively take part in the flow. This is proved by the fact that we can get an expression similar to [START_REF] Sahimi | Non-linear transport processes in disordered media[END_REF] in the case of a series of parallel conduits of different sizes with a maximum diameter equal to D .

The above expression can be extrapolated to any kind of porous medium, under the condition that one is able to identify (i) the appropriate value for the maximum width of the widest path, which will govern the value of the coefficient  , and (ii) the type of pore size distribution and structure, which will govern the value of  .

Figure 1 :

 1 Figure 1: Typical flow curve of a Carbopol solution (squares) and an emulsion (circles). The continuous lines correspond to the HB model fitted to data with the following parameters: (Carbopol) Pa c

Figure 2 :

 2 Figure 2: Functional scheme of the complete set up used for injection (a) and simplified cross-section of the confinement cell (b).

Figure 3 :

 3 Figure 3: Laser granulometry data for glass beads with three different mean diameters ( D ): 0.26 mm (triangles), 0.6 mm (circles), 1.1 mm (squares).

Figure 4 :

 4 Figure 4: Typical set of data for a complete test of injection (bead size 0.26mm, (Carbopol) Pa c 78  

Figure 5 :

 5 Figure 5: Typical interpretation of the mass and pressure data in terms of pressure gradient and flow rate for an injection test of the type shown in Figure 4 (filled squares), and for the same test with the phantom (open squares). The continuous line is a HB model fitted to phantom data (see text).

Figure 6 :

 6 Figure 6: Glycerol injection tests: apparent flow curve for the phantom (after subtracting the hydrostatic pressure drop) (open circles) and effective flow curves for two porous media: mm 1 . 1  D (filled squares), and mm 2  D(filled circles). The continuous lines are linear fits of these data.

Figure 7 :

 7 Figure 7: Effective flow curves for injection tests with a Carbopol solution ( Pa c 78  

Figure 8 :

 8 Figure 8: Effective flow curves for injection tests through a porous medium ( mm 2 D) with

1 VFigure 9 :

 19 Figure 9: Effective flow curves for injection tests through a porous medium ( mm 2  D ) with two Carbopol solutions: (open circles) Pa 55  c

Figure 10 :

 10 Figure 10: Parameters (critical pressure drop (a) and consistency (b)) of the Herschel-Bulkley model (equation (4)) fitted to data of injection tests as a function of the appropriate expressions deduced from dimensional analysis and involving the corresponding rheological parameters of the fluids for Carbopol solutions (squares) and emulsions (stars) and the porous medium characteristics. The continuous lines correspond to straight lines of slope  equal to 12 for the Carbopol and 5.5 for the emulsion in (a), and of slope 85 

Figure 11 :

 11 Figure 11: Different models for representing the pressure drop vs velocity for the steady flow of a yield stress fluid ( 36 . 0  n ,

Figure 12 :

 12 Figure 12: Ratio of the two coefficients of equation (12) as a function of the minimum to maximum diameter ratio for the flow through a series of conduits of different diameters (see text).

Figure 13 :

 13 Figure 13: Mass evolution as a function of time for the flow of Carbopol solution ( c  =65Pa) (open symbols) and emulsion (

that since the Herschel-Bulkley model contains three adjustable parameters one can find slightly different values depending on the fitting procedure used. As in our case we observed that the resulting value for n ranged from 0.34 to 0.38, we decided to fix its value to 0.36 for all materials and adjusted the two other parameters, which makes it possible to represent the data as accurately as through various other procedures.

										8				89	Pa	,	3	.	6		k		49	Pa.s	n	and	n		36 . 0	. For the emulsions we
													c
	found:	54				71	Pa	,	12	7 .		k		13	4 .	Pa.s	n		and	n		36 . 0	. Note
				c