Burst dynamics, up-scaling and dissipation of slow drainage in porous media

We present a theoretical and experimental investigation of drainage in porous media. The study is limited to stabilized fluid fronts at moderate injection rates, but it takes into account capillary, viscous, and gravitational forces. In this theory the work applied on the system, the energy dissipation, the final saturation and the width of the stabilized fluid front can all be calculated if we know the dimensionless fluctuation number, the wetting properties, the surface tension between the fluids, the fractal dimensions of the invasion front and the invading structure, and the exponent describing the divergence of the correlation length in percolation. This theoretical description explains how the Haines jumps' local activity and dissipation relate to dissipation on larger scales.

In most practical applications of porous medium physics, the typical length scales where our interest lies is substantially larger than the scale where the relevant physics is taking place. Oil reservoirs or water aquifers are in the range of kilometers while the typical pore sizes are commonly in the micrometer scale, about nine orders * maloy@fys.uio.no of magnitude smaller. How do we deduce the flow behavior at large scales from small-scale physics? The usual way of solving these problems is a top-down approach using Darcy s law [54] on a mesoscopic level. However this approach does not take into account local fluctuations, like capillary or viscous fluctuations, which are averaged out. In this manuscript we take an alternative bottomup approach where we emphasize on the pore-level capillary fluctuations and compare those fluctuations with the characteristic forces which are set up by the external fields on the whole system. Examples of such forces are gravitational or viscous fields. We will also limit our discussion to the stable drainage regime, which occurs when a nonwetting fluid displaces a wetting fluid, and when the viscous and/or gravitational forces stabilize the displacement front between the two fluids. This approach was introduced in the late 1970s [55,56] and its theoretical development benefited greatly from invasion percolation models [57]. On this subject, other relevant experimental, numerical, and theoretical articles have since been published [26,27,[58][59][60][61]. The dimensionless fluctuation number F , introduced in [30,61], quantifies the ratio between the viscous and/or gravitational field and the capillary pressure fluctuations. The characteristic length scale η which describes the width of the invasion front, and which depends on F , is of central importance to calculate the saturation behind the front, which in turn gives a measure of the sweeping efficiency of a given drainage process, a quantity of great interest in a number of applications. The local structures will be fractal on length scales smaller than η with a crossover to a homogeneous behaviour on length scales larger than η. Knowing the scaling of the fluid structures up to the length scale η allows one to have full control of the energy balance of the problem and calculate large-scale quantities such as the final saturation behind the front, the dissipation or the entropy production, and the work required to move the front forward.
The dynamics of slow displacement in porous media has been observed to occur in an intermittent manner by so called Haines jumps [59,[62][63][64][65][66][67][68][69][70][71][72]. The invasion percolation model [55][56][57] describes well the structure of slow displacement in porous media [22], but it does not describe the dynamics realistically because the invasion is limited to one pore at a time. However it is possible to introduce a realistic interpretation of time in a modified invasion percolation model [63] by introducing a constant κ relating the volume change at the interface with a corresponding change in pressure. Ref. [63,64] found both in simulations and experiments that the pressure fluctuations of the Haines jumps follow a power law distribution with an exponential cutoff function.
In the present work we use a new approach to calculate the dissipation and the energy of the surfaces left behind the front by considering the elastic energy (surface energy) released by the Haines Jumps when the front is in a steady state regime. In between jumps all the external work goes into building up the elastic (surface) energy of the fluid front. When the invading fluid front is in a steady state, its average length is constant, and the work applied to the system (for example by an external pump) must be equal to the elastic energy released by the invasion front (green line in Fig. 1). As a result, using the distribution function of the capillary pressure fluctuations and performing a bottom-up approach to integrate up the elastic energy released by the bursts, we can check the consistency of our theory as well as the distribution of dissipative events. We found as expected that the total elastic energy released by the bursts is equal to the work W = p ∆V , where p is the average pressure across the model and ∆V is the volume change of the invading fluid corresponding to the interface motion. This result provides an important consistency check for the analysis. The dissipation can then be calculated by subtracting the generated surface energy from the total applied work. The surface energy is directly computed by the scaling of the invasion structure, which is given by the invasion structure's fractal dimension, also measured in our work. As a result of the theory and quasi-two-dimensional experiments, we discovered that the work W , dissipation Φ, and the energy of the surfaces left behind the front E s are all proportional to the number of pores S invaded. Up to the characteristic length scale η, the spatial scaling of S is fractal and for larger length scales the scaling becomes proportional to the system size instead. We then show that we can explicitly calculate W , E s , and Φ if we know the surface tension of the fluids and the wetting properties from the experiments. We further discovered an important analytical result: that the ratios E s /W and Φ/W are both independent of system size.

II. EXPERIMENTAL TECHNIQUE
The majority of the experimental results presented in this paper are taken from data published in the past in our group. The system is composed of a modified Hele-Shaw cell [73], filled with a monolayer of glass beads with diameters in the range [1.0mm, 1.2mm]. The glass beads are randomly distributed in the cell gap and the voids in between the beads form the porous network. This quasi-two-dimensional geometry allows for the direct visualization of the fluid phases, by means of regular optical imaging using a digital camera. For about the model construction, see for instance Refs. [19] and [74].
In this paper we also present results produced using an alternative stereolithography 3D printing technique to create the porous structures. We have employed a Formlabs Form 3 printer to produce models in a transparent plastic material (Clear Resin FLGPCL04). This technique allows us to control the geometry of the porous network and in particular to tune its porosity. In the experiments presented here we have made quasi-twodimensional models where cylinders are distributed with a Random Sequential Adsorption algorithm [75] where we can set the minimum distance between the cylinders (see Fig. 1), typically chosen as 0.3mm. The cylinders height and diameter were both chosen to be 1mm. The spatial resolution of these models is about 0.09mm and the models are constructed to optimize the visualization of the pores, which are seen from a top-down view. The 3D printed porous model is placed between two thick Plexiglass plates which are clamped around the edges using screws to give robustness to the setup and ensure the quasi-two-dimensional geometry. A flexible PVC film is placed between the porous model and the top Plexiglass plate. This film plays an important role: due to its flexibility, when the screws around the model are fastened the PVC film gets in contact with the top of all cylinders, thus ensuring the appropriate sealing of the model. The film has similar wetting properties as the 3D printing material used in the construction of the model. Fig. 1 shows a typical snapshot of an experiment and a diagram of the setup. We also define in the upper part of the figure the model's width w, length L, the invading front (green line), its average position h and width 2η.
The porous network is initially fully saturated with a wetting viscous liquid composed of a mixture of glycerol (80% in weight) and water (20% in weight). The kinematic viscosity, density and surface tension (with respect to air) are ν = 4.25 · 10 −5 m 2 /s, ρ = 1.205 g/cm 3 and γ = 0.064 N/m. The wetting liquid is dyed with a dark blue colorant (Nigrosin), to aid visualization. Air is used as the nonwetting phase. The contact angle measured inside the wetting phase is ψ = 70 • . During an experiment, the liquid phase is withdrawn with a syringe pump (Harvard Apparatus) at a constant flow rate, leaving the model from a width-spanning channel at the bottom end of the cell. Air enters the model from the top, through another width-spanning channel that is open to the at-

III. SURFACE ENERGY, DISSIPATION AND BURST DYNAMICS
Consider a single pore identified by the index i in which a nonwetting fluid displaces a wetting fluid as illustrated in Fig. 2. The front moves from a position with a radius of curvature r 0 to one with a radius of curvature r = r 0 −dr. This results in a volume change dv i of the nonwetting fluid, and a corresponding volume change −dv i of the wetting fluid in that pore. Assume that we measure the distances in the normal (radial) direction from each of the points on the front defined with radius r 0 to the front defined by the radius of curvature r = r 0 − dr. The longest of these distances we denote as dx i . The distance dx i to lowest order in dr is given by where depends on the wetting angle ψ (See Fig. 2) and the capillary pressure p 0 given by the Young-Laplace equation Here r is the in-plane, and r 1 the out of plane radius of curvature assumed to be constant in the quasi-two-dimensional experiments. We can then calculate the corresponding increase in capillary pressure dp = γdr Then by using Eqs. (1) and (3) we get The volume change dv i of the invading fluid in the pore can be written as where A i (ψ, p 0 ) is the surface area in pore i which is function of the wetting angle and the capillary pressure p 0 . We then get such that where the capacitive volume gives the fluid volumetric change per unit capillary pressure in pore i [63]. Let U i (x) be the elastic energy (surface energy) of the interface in pore i, and dx i a small displacement of the interface between the two fluids due to an external pump driving the system at a low flow rate such that viscous forces can be neglected. Assume that the system is in mechanical equilibrium. The work performed by the pump will then increase the elastic energy due to this displacement. The elastic energy associated to pore i having interface located at the position x i would then be where p is the capillary pressure across the interface, and U i (0) is the elastic energy when x i = 0. Using Eqs. (7), (9) and (5) we get We have assumed p 0 = 0 without loss of generality since κ i p 2 0 /2 also can be included in U i (0). If instead of considering a single pore we consider the work done on a front having n pores, the total elastic energy U t will be the sum of the contributions from all pores belonging to the invasion front (green line in Fig. 1). We then get where κ = ( where κ i is the capacitive volume of pore i so that κ is the average capacitive volume over the n interface pores. The capillary pressure p 2 at the time t 2 and the capillary pressure p 1 at the time t 1 are related by where q is the imposed volumetric flow rate. Since we have assumed the system to be in mechanical equilibrium, the interface will slowly increase its capillary pressure without any bursts (Haines jumps) [62]. The capillary pressure will then build up from p 1 at t 1 to p 2 at t 2 due to the work W performed by the external pump, see Fig. 3. This work will then increase the total elastic energy (surface energy) of the interface from After some time the interface will however reach a situation where the capillary pressure in one of the porethroats is at the threshold value p t . The interface at that pore will then become unstable, and the invading fluid will move into one or more neighboring pores. At this time, the capillary pressure at the interface of the growing burst will be lower than the capillary pressure at the other parts of the interface. This produces a local velocity field that extends from the growing burst to the other parts of the interface. The interface will then back-contract, beginning with the pores closest to the pore where the burst begins and spreading out through the interface until it reaches equilibrium. Then the capillary pressure is again the same at all pores along the interface. Let us assume that the burst starts at time t 2 at a capillary pressure p 2 = p t and that, after the burst, the system reaches another equilibrium state at capillary pressure p 3 . During the burst the elastic energy of the part of the fluid interface which is back-contracting will be reduced. This energy reduction is equal to and will go to creation of the surface energy of the new burst in addition to viscous dissipation due to the flow that takes place during the burst. Here ∆p = p 2 − p 3 is the pressure drop in the burst and (p 2 + p 3 )/2 is the average pressure during the burst. The surface energy of one single burst in our quasi-twodimensional experiments has two contributions, one from the area of the burst 2γ 2 a 2 s and one from the interface contour aγ 1 bs e . Here a is the average in-plane pore size, b is the height of the cylinders (or the diameter of the beads in the case of the previously published experiments based on a monolayer of glass beads, see Sec. II), s is the number of pores in the burst considered, and as e is the length of the invasion contour of the burst (including trapped clusters). The factor 2 in the term 2γ 2 a 2 s is due to the creation of two new solid-air interfaces after the burst, one at the bottom and one at the top. The need for two different surface energies values, γ 1 and γ 2 , is explained below.
In Fig. 4 a) we show the tracking of the invasion front (green line), and the in-plane invasion contour left behind the front (red line). The fractal dimensions D of the invading fluid phase (air), D e of the invasion contours, and D f of the invasion front are shown in Fig. 4 b), where a box counting technique was employed in the measurement [49]. We have obtained the values D = 1.87 ± 0.10, D e = 1.84 ± 0.10 and D f = 1.45 ± 0.10. The dotted and dashed lines on the left correspond to the exponents -2 and -1 respectively. These are expected as the fractal nature of the invasion structure has a lower bound at the pore-size a. For boxes of size δ < a, we recover the intrinsic 2 dimensional nature for the invading structure mass (dotted line) and the 1 dimensional nature of the perimeters of the front and internal contours (dashed lines).
When a fluid film is left behind the invasion front the surface tension γ 1 = γ 2 = γ which is the surface tension between the two fluids. However when no film is left behind the invasion front, γ 2 = γ ns − γ ws , where γ ns is the interface tension between the nonwetting fluid and the solid and γ ws is the interface tension between the wetting fluid and the solid. Using the Young equation [10] we can also write the previous relation as γ 2 = γ cos(ψ), where ψ is the static contact angle at the liquid-solid-air triple line. In the situation in which no wetting fluid film is left behind the invasion, the in-plane interface contour (red lines in Fig. 4) will be partly formed by nonwettingwetting segments and partly formed by nonwetting-solid phase segments. This division is exemplified in Fig. 5, where we have split the contours shown in red in Fig. 4 into orange and green segments, where orange denotes the interface between the nonwetting and wetting phases and green the interface between nonwetting and solid phases. Let us define the ratio = A ns /A t , where A ns is the interface's surface area between the nonwetting phase and the solid phase, and A t is the total surface area of the interface. In Fig. 5, corresponds to the ratio between the total length of green lines to the total length of both green and orange lines added. We then get an effective surface tension γ 1 = (1 − )γ + (γ ns − γ ws ) for the total internal interfacial contours. Using again the Young equation [10] we can rewrite this expression as We will in the following consider a situation where we have a long channel containing the porous medium with the invading fluid entering from one of the short sides of this channel and with the outlet on the other side. We will further assume that we follow the displacement front. The front will have a width of the order of the width w of the channel. Apart from an initial transient we can then consider the front to be in a steady state regime with a constant average length n (measured in number of pores). Since the elastic energy of the growing interface on average must be constant, the work W applied by the pump must equal the elastic energy released by the bursts which is equal to the dissipation φ plus the surface energy needed to create the the new interface E s , such that It is important to note that W , Φ, and E s are average quantities that do not account for fluctuations, and that φ and E s are averaged over a large number of bursts. From [63] we know that the distribution function of the burst sizes G(s) is directly linked to the distribution of the capillary pressure drops F(∆p) during the bursts since wherev is the average single pore volume averaged over the invading structure, and We found [59,63,64,71], both in simulations and experiments that the burst size distribution G(s) is given by the scaling relation where the cutoff size of the bursts s * is [63] where N (p c ) is the value of the normalized capillary pressure threshold distribution taken at the critical percolation threshold pressure p c , and ν = 4/3 is the exponent describing the divergence of the correlation length in percolation [76]. Since the capillary pressure threshold distribution N (p t ) is normalized, the width of the distribution σ ≈ 1/N (p c ). Therefore the cutoff size of the bursts will increase with σ according to Eq. (20) as s * ∝ σ 0.78 . Martys et al. [77] derived the analytical form where D and D f are, respectively, the fractal dimensions of the growing cluster and its front (seen as the green line Fig. 1). Using the literature values D = 1.82 [22,57] and D f = 4/3 [23,27,58], we obtain τ = 1.32, consistent with our experimental data shown in Fig. 6. We therefore get from Eqs. (17), (18), (19) and (20) where Eqs. (20) and (17) give the cutoff pressure ∆p * Fig. 7 shows the scaling function H(∆p/∆p * ) for various values of κ from invasion percolation simulations [16] which confirms (22) and (23). We will in the following use F(∆p) to calculate the work. Assume that the viscosity is low such that the total time of the bursts is short, and can be neglected The dependence of the crossover function H(∆p/∆p * ) of the pressure jump distribution in a modified invasion percolation simulation [63]. For each value of κ, the points represent averages over five independent simulations on 200 × 1500 lattices. The scaling exponent 0.30 ≈ 1/(1 + νD). Data is taken from [63].
compared to the total time considered. The work can then be calculated by averaging Eq. (15) where M is the total number of bursts during the time considered. Here we have assumed that the average pressure in each burst (p 2 + p 3 )/2 in Eq. (15) is independent of the capillary pressure drop of that burst and therefore can be averaged independently of ∆p. The value of that average is p . This assumption can be readily verified experimentally if we look at the typical temporal signature of the capillary pressure in one experiment, see where I is the integral Then by using Eq. (23) we obtain Assume that we consider the time t w the front needs to move a length w corresponding to the width of the model. Within this time the invading fluid has invaded an area a 2 S of S pores. The total number of bursts considered M can be calculated as S/ s , where s is the average size of the bursts where I is the integral in (26). We then get the following simple expression for the work This expression is the well-known relation of the work expressed as a pressure times a volume change p ∆V , where ∆V = Sv = Sa 2 b = qt w is the volume injected by the pump during the time t w . We can then calculate the total dissipation within t w by using equations (16) and (29) (30) where C is a constant. Here the second term corresponds to the contribution to the surface energy from the top and the bottom interfaces (see discussion after Eq. (15)) and the third term to the contribution from the side walls. The D −1 term in the exponent is due to the cut between the fractal invasion structure and the sidewall using one of Mandelbrot's rules of thumb [48,49]. The last term corresponds to the surface energy connected to the contours of the trapped clusters S e seen in red in Fig. 4. For large S the third term can be neglected and The length of the front n, seen as the green line in Fig. 4, is given by where D f is the fractal dimension of the front (not included trapping) and D e is the fractal dimension of contours of the trapped clusters. From Fig. 4 we see that the measured fractal dimensions for D e and D seem very close, effectively indistinguishable from one another considered the error bars. This result implies that S ∝ S e . In the experiment shown in Fig. 4 we have measured the proportionality relation to be S e = 1.32S. We have also measured the prefactor linking S and the window size w/a (see relation (32)), obtaining the expression S = 0.72(w/a) D . By plugging these two equations into Eq. (31) for the dissipation φ, we gain access to the full energy balance of the system, the resulting plot is shown in Fig. 9. In this figure we see how the total applied work W (thick black line) is split into the dissipated energy φ (blue) and the surface energy associated to creation of the top and bottom interfaces E st = γ 2 a 2 2S (green) and the lateral interface E sl = γ 1 abS e (yellow). The dissipation is calculated directly from Eq. (31). The total applied work, the surface energy associated to the top and bottom interfaces and that associated to the lateral interfaces are given respectively by the first, second and third terms in Eq. (31). Because all terms in (31) and (29) are proportional to S the ratios φ/W , E sl /W and E st /W are independent of w/a for large systems.
In the computation of the energy balance in Fig. 9 using Eqs. (31) and (32), we used = 0.71 measured from the experimental images (see Fig. 5), γ = 0.064N/m, ψ = 70 • , γ 1 = (1 − )γ + γ cos(ψ) = 0.034N/m, p = 210P a and D = 1.87, see Fig. 4. As previously stated, our assumptions for γ 1 and γ 2 correspond to the scenario in which no film of the defending wetting fluid is left behind covering the solid surfaces after the invasion by the nonwetting phase.
If instead of the quasi-two-dimensional system considered thus far we had a three-dimensional long square channel with a cross sectional area of w 2 , Eq. (29) for the work would remain the same, while Eq. (30) for the dissipation would be modified to where C is a constant. When the system gets large the last term can be neglected compared to the two first terms such that If the invaded surface area S e a 2 times the characteristic pore size a is proportional to the invaded volume Sa 3 , the ratio φ/W and E sl /W will be independent of w/a.

IV. STABILIZING FIELDS AND CROSSOVER LENGTHS
Let's now take into account a gravitational field and consider an experiment where a low density fluid with density ρ 1 and viscosity µ 1 = 0 invades another fluid with a higher density ρ 2 and viscosity µ 2 from above. (See Fig. 1). The mapping between the occupation probability in percolation theory f and the capillary pressure p is given by [76,78] We therefore have where the critical occupation probability is f c and the critical percolation pressure is p c . As before, N (p t ) is the normalized capillary pressure threshold distribution. Now, by Taylor expanding N (p t ) around p c to the lowest order in p − p c in Eq. (36) we get It is reasonable to assume that the viscous pressure drop in the displaced fluid will depend linearly on the length scale since there are no trapped invading fluid clusters in the displaced fluid. Let us consider the capillary pressure at a height x. Since the gravitational field also depends linearly on the length scale x we have where x 0 corresponds to a the height along the front where the capillary pressure is at the percolation threshold. Here k is the permeability felt by the displaced fluid, A the cross-section area of the porous medium, and ∆ρ = ρ 2 − ρ 1 the density difference between the fluids. The fluctuation number F is a dimensionless number [61] which characterizes the gravitational and viscous fields and the capillary pressure fluctuations. Note that the width of the capillary pressure threshold distribution is σ ≈ 1/N (p c ), as the distribution is normalized. We will then use an assumption, first introduced by Sapoval [27,57,79,80], that the correlation length in percolation will scale as the width of the front η = (x − x 0 ) ∝ ξ, where ξ is the percolation correlation length. This is based on the observation that the largest trapped cluster of wetting fluid is limited by the front width. From percolation theory, the correlation length ξ is given by [76] (ξ/a) Using Sapoval's argument, and inserting Eq. (40) into Eq. (38) gives [26,61] Therefore, when F > 0, the front will be stabilized with a characteristic length scale η [26,27,61] . This scaling behaviour is shown in Fig. 10. The blue dots show experimental data taken from [61] where the authors analyzed the scaling of the invasion front in drainage by keeping the Bond number Bo = ∆ρga 2 /γ = 0.154 fixed but changing the capillary number Ca by varying the flow rate q (see Fig. 10 in Ref. [61] and experiments 1, 2, 3, 4, 5, 7, and 9 in Table I in the same reference). Let us first assume that the gravitational effect is large enough such that it is the characteristic length scale η that sets the width of the front η < w. Consider again a time scale t w corresponding to the time the front needs to move a length scale w. Dividing space into cells of size η 2 , below which the fluid distribution is fractal, the number of pores S and S e are given by and the work can be written as where h is the average position of the front from the lower outlet. The first term in (43) is the negative work due to the hydrostatic pressure, the second term the work due to the viscous pressure in the displaced fluid, and the last term is the work needed to continuously build up the interface energy. The work must be equal to the dissipation Φ plus the surface energy E s minus the change in gravitational potential energy ∆U p (corresponding to the first term in Eq. (43)) We then get the dissipation which for large S can be approximated as Note that in the limit of very slow processes, the flow rate q becomes negligible so the first term in this Eqs. (45) and (46) can be ignored, thus recovering Eqs. (30) and (31).
In the scenario in which η > w, the relations for S and S e will change to but we will still have the same expressions for W and φ. Again if D e = D, the ratios φ/W , E sl /W and E st /W will be independent of w/a for large systems.

V. SATURATION BEHIND THE INVASION FRONT
The upscaling problem of calculating the large scale saturation involves identifying the cross-over length scale where the fluid distribution is no longer fractal. Let L be the length of the porous model and assume L > w. The final saturation behind the front of the invading fluid and its dependence on the pressure across the model has been studied in Ref. [81]. In this study, we considered the volume of the invaded fluid in boxes with size corresponding to the width of the invading front. On length scales below this size, the structure of the invading fluid is fractal, while on length scales larger than this size, the structure is homogeneous. For sufficiently large F , when η is the characteristic length scale of the front, the saturation S F nw of the nonwetting fluid behind the front becomes [81] S F nw ∝ where d is the spatial dimension (2 or 3). Hence using Eq. (41) However, when the fluctuation number F is sufficient small, η > w, and the width of the model w will be the characteristic length scale in the problem. Then Fig. 11 shows a two-dimensional invasion percolation simulation with a gravitational field together with drainage experiments performed by Ayaz et al. [81]. The red dash-dotted line confirms the predictions of the theoretical scaling in Eq. (49). We make the important remark here that in the equations considered above we have neglected the boundary effects at the inlet and outlet since L/w 1. These boundary effects are however important when the characteristic length scale of the front w is of the same order as the length L of the system [74].
In section IV we considered an invading fluid with a negligible viscosity µ 1 = 0. Let us now instead consider the scenario in which this fluid has a non-negligible viscosity µ 1 . In that case the fluctuation number will be where k 1 is an effective permeability that depends on the characteristic crossover length scale of the system which is either η or w. The expression (51) is meaningful only if the pressure drop in the invading fluid depends linearly on the length scale along the average flow direction. Network simulations in two dimensions however, show a linear change in the capillary pressure due to viscous flow [60], [59]. This is due to the loop-less strands in the invading nonwetting fluid. It is of great interest to know if this is also valid in three dimensions. The structure left behind the front is a fractal structure with clusters of all length scales up to the characteristic length scale η or w [27,58]. Therefore if η < w, the permeability k 1 (F ) will be a function of F . If we know this functional dependence we can solve the equation with respect to F and thereby find the characteristic length scale η from We then also know the saturation through Eq. (48). A goal for future experiments/simulations/theory should therefore be to measure/predict the permeability dependence of the invading fluid as function of the characteristic length scale η.

VI. CONCLUSION
In this paper, we discussed the importance of capillary fluctuations in porous media, as well as the characteristic length scales set by the competition between capillary fluctuations and external fields, gravitational or viscous. We focused on fluid fronts that are stabilized by gravitational and viscous fields. The fluctuation number F , which describes the scaling of the front width η, was introduced. When considering viscous and gravitational fields, the theory describes well the scaling of the width of the fluid front and the final saturation of the fluid left behind the invasion front observed in experiments. On a length scale smaller than η, the structure within the front is generally fractal, while on a length scale larger than η, it is homogeneous. The characteristic length scale η is thus of primary importance in defining a relevant Representative Elementary Volume (REV) for an average Darcy description of the two-phase flow problem.
We also discussed the energy dissipation caused by out-of-equilibrium Haines jumps and calculated the total elastic energy released by those jumps using the pressure drops connected to the jumps' known scaling behavior. When the invasion front is in a steady state regime, we find that the elastic energy released by the jumps corresponds to the work applied to the system, as expected. However, this description also explains the local activity and dissipation caused by local Haines jumps, as well as how it is related to dissipation on a larger scale.
We calculated the energy dissipation and saturation left behind the invasion front in addition to the work done by the pump to push the fluid front forward. These quantities are affected by the generalized fluctuation number, as well as the fractal dimension of the invasion structure, the fractal dimension of the contour of the trapped clusters, and the exponent ν, which describes the divergence of the correlation length in percolation as it approaches the critical percolation threshold f c . The work W , dis-sipation Φ, and surface energy E s were discovered to be proportional to the number of pores invaded S, with fractal scaling on length scales smaller than the characteristic length scale, η or w, but scaling with the spatial dimension on larger length scales. Hence, the ratio of these energies stays constant as the system size increases.
The theory described here is a bottom-up approach that allows the problem to be scaled up from small to large scales. The theoretical results were compared to and found to be consistent with quasi-two-dimensional experiments. Because of the small temperature increase involved in heat dissipation from the Haines jumps (mK or less), it is very challenging to measure these temperature fluctuations. Such a measurement remains to be done. An experimental verification of the theoretical scaling of the invasion front and the saturation behind the front for three-dimensional porous media is also of great interest for future work [81]. Experiments involving changing the gravitational effects in a geophysical centrifuge could be one way to accomplish this. The porous medium in a real rock is typically not homogeneous, in contrast to the porous media considered here. As a result, it is of great interest to extend the theory to include also the case of inhomogeneous porous media, for example with a porosity gradient.