Table Of Content

Influence of pore-scale disorder on viscous fingering during drainage Renaud Toussaint,1 Grunde Løvoll,1 Yves Méheust,2 Knut Jørgen M˚aløy,1 and Jean Schmittbuhl3 1Department of Physics, University of Oslo, POBox 1048 Blindern, N-0316 Oslo, Norway 2Department of Physics, NTNU Trondheim, N-7491 Trondheim, Norway 3Laboratoire de Géologie, Ećole Normale Supérieure, 24 rue Lhomond, F-75231 Paris, France (Dated: February 2, 2008) We study viscous fingering during drainage experiments in linear Hele-Shaw cells filled with a 5 random porous medium. The central zone of the cell is found to be statistically more occupied 0 0 than the average, and to have a lateral width of 40% of the system width, irrespectively of the 2 capillary number Ca. A crossover length wf ∝ Ca−1 separates lower scales where the invader’s fractaldimension D≃1.83 isidenticaltocapillary fingering,andlarger scaleswherethedimension n isfoundtobeD≃1.53. Thelateral widthandthelarge scale dimensionarelower thantheresults a for Diffusion Limited Aggregation, but can be explained in terms of Dielectric Breakdown Model. J Indeed,weshowthatwhenaveragingoverthequencheddisorderincapillarythresholds,aneffective 4 law v∝(∇P)2 relates theaverage interface growth rate and thelocal pressure gradient. ] t PACSnumbers: 47.20.Gv,47.53.+n,47.54.+r,47.55.-t,47.55.Mh,68.05.-n,68.05.Cf,81.05.Rm. f o s . Viscous fingering instabilities in immiscible two-fluid cells belong to the family of Laplacian Growth Models, at flows in porous materials have been intensely studied i.e. obey the Laplacian growth equation ∇2P = 0, with m over the past 50 years [1], both because of their impor- an interfacial growth rate v ∝−∇P, where P is the dif- - tant role in oil recovery processes, and as a paradigm of fusing field, i.e. the probability density of random walk- d simple pattern forming system. Their dynamics is con- ers in DLA, or the pressure in viscous fingering. Despite n trolled by the interplay between viscous, capillary and differences as respectively a stochastic and deterministic o c gravity forces. The ratio of viscous forces to the cap- growth, and boundary conditions as respectively P = 0 [ illary ones at pore scale is quantified by the capillary orP =−γ/rwithr theinterfacialcurvature,itisadmit- number Ca = µv a2/(γκ), where a is the characteristic ted that these processes belong to the same universality 1 f v pore size, vf is the filtration velocity, γ the interfacial class [2, 3, 4]. In radial geometry these processes lead 7 tension, κ the permeability of the cell, and µ the viscos- to fractal structures of dimensions D = 1.70±0.03 [5], 4 ityofthedisplacedfluid,supposedheremuchlargerthan 1.713±0.0003 [6], and 1.7 [4] respectively in viscous fin- 0 the viscosity of the invading one. gering in empty Hele-Shaw cells, DLA, and numerical 1 0 There is a strong analogybetween viscous fingering in solutions of deterministic Laplacian growth. The two 5 porousmediaandDiffusionLimitedAggregation(DLA), numerical models have been reexplored recently using 0 aswasfirstpointedoutbyPaterson[2]. Indeed,bothpro- stochastic conformal mapping theory [3, 4, 6]. However, t/ cesses of DLA and viscous fingering in empty Hele-Shaw in Hele-Shaw cells filled with disordered porous materi- a alssimilartotheoneusedhere,alowerfractaldimension m D =1.58±0.09 has been measured [7]. - d In straight channels, DLA gives rise to fractal struc- n tures of dimension 1.71, occupying on average a lateral o fraction λ = 0.62 of the system width W [8]. Viscous c fingering in empty Hele-Shaw cell converge towards the : v SaffmanTaylor(ST) solution[9],with auniformly prop- i X agating fingerlike interface covering a fraction λ = 0.5 of the system width at large capillary numbers [9, 10], r a selected by the interfacial tension [11]. In the system we study, the cell is filled with a disor- deredporousmedium,andthenon-wettinginvaderoflow viscosityshowsabranchedstructurethatdependsonCa (Fig.1). Wewillshowtheoreticallythatifindeedathigh capillary number, the process is well described by DLA FIG. 1: Invasion clusters on thresholded images at capil- as often suggested [2], there is also at intermediate Ca a lary numbers Ca = 0.06 (a) and 0.22 (b,c), for W/a = 210 (a,b)and 110 (c), with displayed system lateral boundaries. regime where the flow in random porous media is better The superimposed gray map shows the occupancy probability described by another Laplacian model, namely a Dielec- π(x,z)oftheinvader,inamovingreferenceframeattachedto tricBreakdownModel(DBM)withη =2–theinterfacial the most advanced invasion tip and to the lateral boundaries. growth rate is v ∝(∇P)η in DBM, η =1 corresponding 2 haveshownthat the invasionprocessis stationary,upto 11 fluctuationsarisingfromthedisorderinporegeometries. To extract the underlying average stationary behavior, 00..88 Ca all quantities are then analyzed in the reference frame n(z)/n 0000....4466 W/a = 110 W/a = 215 W/a = 430 (xx=,z)Wa,ttaancdhetdo tthoetfhoerelmatoesrtalpbroopuangdaatriniegstaipt xat=z =0 a0n,dz Ca = 0.03 00..22 ΦST(z, )λ=0.4 CCCaaa === 000...012622 pagoaintteisngatagaarinosutgthhlye flcoownstdainretctsipoened(thvitsipti[p14i]n)d.eeAdvperraogpe- quantities at any position (x,z) of the tip related frame, 00 00 0.5 11 1.5 22 2z/W aredefinedusingallstagesandpointsoftheinvasionpro- cess,excluding regionscloserthanW/2 fromthe inlet or FIG. 2: Scaled invader’s density n(z)/nCa as function of outlet, to avoid finite size effects. 2z/W, distance to the most advanced tip scaled by the sys- The average occupancy map π(x,z) is defined as [14]: tem’shalfwidth,forthreesystemsizesandfourcapillarynum- for each time (or each picture), we assign the value 1 to bers. The lines are the experimentally determined cumulative growthprobability,andthetheoreticalSaffmanTaylorsolution the coordinate (x,z) if air is present there, 0 otherwise. for a single finger occupying a lateral fraction lambda = 0.4 π(x,z) obtained as the time average of such occupancy of the system. function, is displayed as graymap in Fig. 1. Next we compute the average number of occupied poresperunitlengthatadistancez behindthetip,n(z), to DLA –. Our conclusion is supported by both experi- whichisrelatedtoπ asn(z)=(1/a2) W π(x,z)dx. We mentalevidencebasedonacharacterizationofthelarge- 0 showinFig.2adatacollapsefordifferRentcapillarynum- and small-scale geometry of the invader, and theoretical bers Ca and system widths W, n(z)/n = Φ(2z/W), arguments based on averaging the capillary forces con- Ca wheren =(W/a2)v /v [14]. Φisafunctionincreas- tribution to the interface growthrate overthe pore scale Ca f tip ing from 0 at z = 0 towards 1 at z = +∞, as granted randomness. Thevanishingcapillarynumberlimitofour by conservation of the displaced fluid for a statistically system corresponds to capillary fingering, where pores stationary process [14]. Φ evaluated in Fig. 2 is a cu- are invaded one by one, forming a propagating front fill- mulative growth probability defined in Ref. [14], and is ing the whole system, leaving behind isolated clusters of obtained as an average over all experiments and sizes. viscousfluid. Theinvader’sfractaldimensionisD =1.83 [12], theoretically explained by the invasion percolation We also characterize the lateral structure of the in- model[13]. Atmoderatebutfinitecapillarynumbersun- vader in the frozen zone, z > W, where less than der interesthere, we will show that small scalesstill cor- 10% of the invasion activity takes place since Φ(2) > respond to invasion percolation, up to a crossover scale 0.9. We define over this zone a distribution ρ(x) = wf/a ∝ Ca−1 characterizing the maximum size of iso- [W/(a2nCa)] π(x,∞) = (vtip/vf) π(x,∞), so that W lated clusters of defending fluid, abovewhich the system ρ(x)dx/W ≃1. Thisquantity,presentedinFig.3(a) 0 geometry can be described by the DBM with η =2. fRoranaverageoverfiveexperimentsatcapillarynumbers We study viscous fingering processes in linear Hele- Ca=0.06 and 0.22 for W/a=215,and over four exper- Shaw cells of thickness a = 1mm filled at 38% with iments with 0.06<Ca<0.22 for W/a =110, is reason- a monolayer of randomly located immobile glass beads ablyindependentofthecapillarynumberandthesystem of diameter a, between which air displaces a solution size,thoughthenoiseislargerathighestspeed. Thefrac- of 90% glycerol - 10% water of much larger viscosity tion λ of the system, occupied by the invader at satura- µ = 0.165 Pa·s, wetting the beads and walls of the cell, tion is evaluated as in [8]: λ = 1/ρmax, or alternatively i.e. in drainage conditions. The interfacial tension and λ = (x+ − x−)/W, where ρ(x+) = ρ(x−) = ρmax/2. the permeability of the cell are respectively γ = 0.064 Both definitions lead to λ ≃ 0.4±0.02 for the capillary N·m−1 and κ = 0.00166± 0.00017 mm2. We investi- numbers probed, as shown in Fig. 3(a), which is signif- gate regimes ranging from capillary to viscous fingering icantly smaller than the off-lattice DLA result λ = 0.62 (0.01 < Ca < 0.5), in cells with impermeable lateral [8, 15]. walls and dimensions W ×L×a, with widths perpen- The 2D occupancy map π(x,z) itself, displayed as dicularly to the flow direction W/a = 110, 215 and 430, graymap in Fig. 1, has a maximum π =ρ v /v , max max f tip and a length L/a = 840. The cell is set horizontally, so alongalineat(x=W/2,z >W). SimilarlytoArneodo’s that gravity is irrelevant. A constant filtration rate of procedure [8], we determine the support of π > π /2, max water-glycerol is ensured by a controlled gravity-driven displayed in Fig. 3(b) and (c), which corresponds to the pump. most often occupied region. Within the noise error, the Picturesoftheflowpatternaretakenfromthetop,and shape of this region resembles the theoretical ST curve treatedto extractthe invading air cluster (with pixels of corresponding to this λ = 0.4 [9] (gray lines in Fig. 3). size0.55a),astheblackclustersinFig.1. Inref.[14],we For such ST finger, this curve would also correspond to 3 4 W/a=215, Ca=0.06 variationsarelowerthancapillarythresholdfluctuations, 3,5 W/a=215, Ca=0.22 W.n(x)/(a.n)Ca2,253 ρmax Wno/ram=a1l1iz0e, d0 .s0q6u<aCrea f<u0n.c2t2ion aefisnntdgertrahinnedgmo,mtohsuttshlirlkeeeaslhdyoinilndgvsta,odweDdhi=pcho1r.ec8so3rc.roeCrsrpoenospvneodrnssdetlotyo,catathpeliallrloagwreyr- ρ (x) = 01,,551 (a0.)5 ρmax tsiiazellsy(sby>twhef)s,ptahteiailnvvaasriioantioanctsivoifty∇iPs d.etAersmsuinmeidngestsheant- b 00 0,2 0,4 0,6 0,8 1 ∇P scales as the imposed ∇P(−∞)∝ Ca, i.e. neglect- x/W b ing the geometryvariationsbetweendifferentspeeds, w f scales as W /∇P ∝ a/Ca, as confirmed by the data t b collapse in Fig. 4(a). w can also be determined experi- (b) (c) f mentallyasacharacteristicbranchwidth,sincecapillary fingeringleavesisolatedclustersoftrappedfluid,whereas FIG. 3: (a) Normalized occupation density ρ(x), with half- the large scale structure is branched. After removing all maximum reached over a width 0.4W. (b,c) Average occupa- trapped clusters, we determine w as the averagelength f tiondensity mapoftheinvaderthresholded athalf-maximum, of intersects of the structure from cuts along x. Indeed, for a system size W/a=215, in the reference frame attached Fig. 4(b) is consistent with a scaling law w /a ∝ Ca−1, to the tip position, at Ca = 0.06 (b) and 0.22 (c), compared f below a saturation at large Ca. to ST curves for λ=0.35 and 0.45. Eventually, we sketch a possible explanation for the 10000 width selection λ = 0.4 and the large scale fractal di- 2 mension D = 1.53±0.02, which are smaller than their 1.8 mber of boxes 100 000C...000a112:497 Slope -1.83 log(w/a)10f111...2461-2 (b) -1.5 C-1a-1 -0.5 cgaoleuLcnattipnelrgapctaiharnetsspmirneaslDsluLsrcAea,lfiereepldsepriemncteitavhbeeillyidtye0f.ev6na2driiaanntgidofln1us.7ild1e.a.dTsNhteoe- Normalized nu 1 00000000........000011223458171976977090 log10(SClao)pe -1.53 c−ibseo∞sutshn)eids=autrnhy−ietcnµovvnefednc/ittκtioir,orenlayanflodocrno∇ngthtPrexo.(plxlreeTd=shsbue0yr,edWtyfihnee)ald·mbxîoiscu=stnhdoe0afnrtyw∇hhPeceo(rpnzerdo=ixˆ-- 0.01 0.361 tion along the invading fluid, i.e. by the capillary pres- (a) sure drop across the meniscus in the pore neck and the 0.01 0.1 1 10 100 (Scale/a).Ca pressure gradient in the invaded fluid. For a given pressure difference at pore scale between the invading air at FIG.4: (a)Massfractal dimensionoftheinvasioncluster at P , and the pressure P in the glycerol-filled pore, we various capillary numbers. (b): Cross-over scale w between 0 1 f capillary and viscous scales, as function of Ca. decompose P0 − P1 = ∆Pv + Pc, where ∆Pv is a viscous pressure drop in the pore neck, and the capillary pressure drop is P = γ/r + 2γ/a, where the in- and c the scaledlongitudinal density of invader Φ(z). Φ(z) de- out-of-plane curvature of the interface are respectively r termined for our experiments (continuous line in Fig. 2) and a/2. As a meniscus progresses between neighbor- and this theoretical ST shape (dashed line) also present ing beads, its curvature goes through a minimum r in m some similarities. Systematic deviations from the math- the pore neck. The meniscus will be able to pass the ematical ST solution at λ=0.4 might nonetheless exist, neck if the pressure drop P − P exceeds the thresh- 0 1 since they have been seen between the DLA envelopes old P = P (r ). For the sake of simplicity, the prob- t c m and the corresponding ST solution at λ=0.62 [15]. ability distribution of the thresholds g(P ) is considered t The mass fractal dimension of the invasion clusters is flat, between P and P , with W = P −P min max t max min analysed by box-counting. N(s) is the number of boxes and g(P ) = θ(P −P )θ(P −P )/W , where θ is t t min max t t ofsizestocovertheinvader. Fig.4displaysanormalized the Heaviside function. In the pure capillary fingering distribution N(s)/N(a/Ca) as function of (s/a)·Ca for limit Ca → 0, the pressure field P is homogeneous in various capillary numbers. By linear regression of this the defending fluid, and a pore is invaded when P −P 0 collapsed log-log data, we find that N(s) ∼ s−1.83±0.01 reachestheminimumthresholdalongtheboundary,close for small scales s < a/Ca, and N(s) ∼ s−1.53±0.02 for to P . At higher capillary number, we want to relate min larger scales. The result can be explained by the fol- the invasion rate to the local capillary threshold, and to lowing approximations. The distribution of pore throats thepressureP intheliquid-filledporenearesttothein- 1 sizesresultsinadistributionofcapillarypressurethresh- terface. If P −P <P , the meniscus adjusts reversibly 0 1 t oldsP ,g(P ),ofcharacteristicwidthW . Considerabox in the pore neck, and the next pore is not invaded. Con- t t t of size w along the cluster boundary in the active zone, versely, if P −P > P , the pore will be invaded, and f 0 1 t such that w ·∇P = W where ∇P is a characteris- most of the invasion time is spent in the thinnest region f b t b tic pressure gradient. At scales s<w , viscous pressure of the pore neck. A characteristic interface velocity can f 4 beevaluatedbytheWashburnequation[16]atthispoint: power-law effective relationship. It would be interesting v ∼ −(2κ/µa)(P −P −P )θ(P −P −P ), where the in future work to explore numerically and experimen- 0 1 t 0 1 t heaviside function results from a zero invasion velocity tally the detailed effect of non-flat capillary threshold if the pore is not invaded. Hypothesizing that only the distributions on the selected fractal dimension, average average growth rate controls the process, independently width occupied in the system, and total displaced mass of the particular realization of random thresholds, the n (Ca) (reported in [14] for the present work), to ex- Ca growth rate averaged over all possible pore neck config- tract the influence of the disorder on the best capillary urations, is number to select in order to maximize the efficiency of the extraction process. −2κ <v >= (P −P −P )θ(P −P −P )g(P )dP Weacknowledgewithpleasurefruitfuldiscussionswith Z µa 0 1 t 0 1 t t t E. G. Flekkøy, A. Lindner, A. Hansen and E. Somfai. = −(κ/µa)θ(P −P −P )× (1) 0 1 min This work was supported by the CNRS/NFR PICS pro- {[(P −P −P )2/W ]θ[P −(P −P )]+ gram, and the NFR Petromax program. 0 1 min t max 0 1 2[P −P −(P +P )/2]θ[P −P −P ]}. 0 1 min max 0 1 max Atmoderatecapillarynumbers,suchasP −P <P , 0 1 max if we assume that the capillary pressure drop is around [1] D. Bensimon, L. P. Kadanoff, S. Liang, B. I. Shraiman, P =P whentheinvasionmeniscusisattheentrance c min and C. Tang, Rev.Mod. Phys. 58, 977 (1986). of the pore neck, we note that (P − P − P )/a = 0 min 1 [2] L. Paterson, Phys. Rev.Lett. 52, 1621 (1984). ∆Pv/a ∼ ∇P/2, and Eq. (1) implies that the growth [3] M. G. Stepanov and L. S. Levitov, Phys. Rev. E 63, rate goes as < v >= −aκ/(4µWt)(∇P)2. This effective 061102 (2001). quadratic relationship between the average growth rate [4] A.LevermannandI.Procaccia,Phys.Rev.E69,031401 and the local pressure gradient arises from the distribu- (2004). tion of capillary thresholds, and means that such inva- [5] E. Sharon, M. G. Moore, W. D. McCormick, and H. L. Swinney,Phys. Rev.Lett. 91, 205504 (2003). sion process should be in the universality class of DBM [6] B. Davidovitch, A. Levermann, and I. Procaccia, Phys. with η = 2, rather than DLA (DBM, η = 1). Indeed, Rev. E62, R5919 (2000). in DBM simulations in linear channels, λ is a decreas- [7] E. L. Hinrichsen, K. J. M˚aløy, J. Feder, and T. Jøssang, ing function of η (as in related deterministic problems, J. Phys. A 22, L271 (1989); K. J. M˚aløy, J. Feder, and as viscous fingering in shear-thinning fluids [17], or η- T. Jøssang, Physical Review Letters 55, 2688 (1985). model [18]), and Somfai et al. [15] report λ ≃ 0.62 and [8] A. Arneodo, J. Elezgaray, M. Tabard, and F. Tallet, Phys. Rev.E 53, 6200 (1996). 0.5 for respectively η = 1 and 1.5, so that the observed [9] P.G.SaffmanandG.Taylor,Proc.Soc.LondonSer. A λ = 0.4 is consistent with η = 2. The fractal dimension 245, 312 (1958). ofDBMisalsoadecreasingfunctionofη,andη =2cor- [10] P.Tabeling,G.Zocchi,andA.Libchaber,J.FluidMech. respondsD =1.4±0.1[19],whichisclosetotheobserved 177, 67 (1987). D =1.53±0.02 in our experiments. [11] J. W. McLean and P. G. Saffman, J. Fluid Mech. 102, Note that at high capillary numbers such that locally 445 (1981); B. I. Shraiman, Phys. Rev. Lett. 56, 2028 P −P ≫P ,thethresholdfluctuationsarenotfeltby (1986). 0 1 max theinterface,andEq.(1)leadsto<v >=−(κ/µ)∇P, [12] R. Lenormand and C. Zarcone, Physical Review Letters inv 54, 2226 (1985). which would correspond to a classic DLA process. We [13] R.Chandler,J.Koplik,K.Lerman,andJ.F.Willemsen, have checked by numerically solving the Laplace equa- Journal Fluid Mechanics 119, 249 (1982). tion with the experimental clusters as boundaries that [14] G.Løvoll,Y.Méheust,R.Toussaint,J.Schmittbuhl,and allexperimentsperformedherewereatmoderateenough K. J. M˚aløy, Phys. Rev.E 70, 026301 (2004). capillary number to have P − P < P all along [15] E. Somfai, R. C. Ball, J. P. DeVita, and L. M. Sander, 0 1 max the boundary [14], i.e. the quadratic law < v >= Phys. Rev.E 68, 020401(R) (2003). −aκ/(4µW )(∇P)2 is expected to hold. [16] E. W.Washburn,Phys. Rev. 17, 273 (1921). t [17] A. Lindner, D. Bonn, E. Corvera Poiré, M. Ben Amar, EvenatmoderateCa,deviationsfromtheDBMmodel and J. Meunier, J. Fluid Mech. 469, 237 (2002). with η = 2 could be observed for significantly non-flat [18] M. Ben Amar, Phys.Rev.E 51, R3819 (1995). distribution of the capillary thresholds in the random [19] J. Mathiesen and M. H. Jensen, Phys. Rev. Lett. porous medium, for which Eq.(1) would lead to a more 88, 235505 (2002); E. Somfai et al. (2004), cond- complicateddependence ofthe growthrate v on∇P, re- mat/0401384; L.PietroneroandE.Tossati,eds.,Fractals flecting the details of this distribution, and not simply a in Physics (North-Holland, 1986), p. 151.