Immiscible two-phase flow in porous media produces different types of patterns depending on the capillary number Ca and viscosity ratio M. At high Ca, viscous instability of the fluid–fluid interface occurs when the displaced fluid is the more viscous, and leads to viscous fingering, which is believed to exhibit the same growth behavior as the viscously-unstable fingers observed in Hele–Shaw cells by Saffman and Taylor [“The penetration of a fluid into a porous medium or Hele–Shaw cell containing a more viscous liquid,” Proc. R. Soc. London 245, 312 (1958)], or as diffusion-limited aggregates (DLA). In such Laplacian growth processes, the interface velocity depends linearly on the local gradient of the physical field that drives the growth process (for two-phase flow, the pressure field). However, a non-linear power-law dependence between the flow rate and the global pressure drop, reminiscent of what has also been observed for steady-state two-phase flow in porous media, was evidenced experimentally for the growth of viscously-unstable drainage fingers in two-dimensional porous media, 20 years ago. Here, we revisit this flow regime using dynamic pore-network modeling and explore the non-linearity in the growth properties. We characterize the previously unstudied dependencies of the statistical finger width and non-linear growth law's exponent on Ca, and discuss quantitatively, based on theoretical arguments, how disorder in the capillary barriers controls the growth process' non-linearity, and why the flow regime crosses over to Laplacian growth at sufficiently high Ca. In addition, the statistical properties of the fingering patterns are compared to those of Saffman–Taylor fingers, DLA growth patterns, and the results from the aforementioned previous experimental study.

Fingering patterns are one of many unique features of two-phase flow, which are caused by hydrodynamic instabilities between two fluids.1,2 When one fluid displaces another fluid inside a medium, depending on the properties of the two fluids and the medium in which they are flowing, the displacement front may exhibit fingering instead of a stable interface between the two fluids.3–5 A wide variety of fingering patterns are observed for different types of multiphase flow, such as the flow of miscible6–9 and immiscible10,11 fluids in continuum1,12 and porous medium,13–15 reactive transport flow,16 and the flow of frictional fluids in granular materials.17,18 The structural properties of different types of fingers are controlled by the underlying physical forces, such as the viscous, capillary, inertial and frictional forces. Furthermore, depending on the driving condition or the geometry of the system, the displacement front can also undergo a transition from fingering patterns to compact19 or foam-like structures.20 

Our study here deals with non-reactive two-phase flow in a non-deformable medium, where the two fluids are immiscible and separated by interfaces associated with a surface tension. Inside a continuum medium, for example, a Hele–Shaw cell, which consists of two parallel plates separated by a small gap, the displaced fluid having a higher viscosity than the displacing fluid leads to viscous instability of the front and the development of fingers with a constant width, except in the vicinity of the fingertip. The smooth rounded shape of these fingers in a two-dimensional (2D) zx plane can be described by the following parametric equation given by Saffman and Taylor in their pioneering work in 1958:12,
z ( x ) = W ( 1 λ ) 2 π ln [ 1 2 ( 1 + cos 2 π x λ W ) ] ,
(1)
where z is the direction of the overall front propagation, W is the width of the flow cell, and the parameter λ is the ratio of the width of the finger to the width of the channel. Note that this theory does not prescribe λ and that the surface tension between the two fluids does not enter this equation, as it was not taken into account in its derivation. However, Saffman and Taylor experimentally measured values of λ close to 0.5 over a wide range of flow rates, and it was later found that the finger width was selected by surface tension.21 When the flow takes place inside a porous medium, on the other hand, for example, a Hele–Shaw cell randomly filled with glass beads, the fingers exhibit a fractal structure.22 Properties of these fractal fingers are controlled by the competition between viscous and capillary forces, the ratio of which is called the capillary number, and the viscosity contrast of the two fluids.23–25 For slow injection, the process is controlled by the disorder in the capillary forces at the pores, and the displacement front generates capillary fingers with a fractal dimension similar to invasion percolation clusters.26 For fast injection of a lower-viscous fluid into a higher-viscous fluid inside a porous medium,27 the growth of the front is governed by viscous instabilities, and it produces fractal viscous fingers with a lower fractal dimension.15,28 It was pointed out that the statistical properties of viscous fingers in porous media are analogous to diffusion-limited aggregation (DLA),29 which is a process of aggregation of matter limited by the diffusion of random walkers arriving from a far distance with a steady flux.30 Both Saffman–Taylor viscous fingering and DLA follow Laplacian growth:
2 p = 0 ,
(2)
where p is the probability density of the random walker for DLA or the pressure drop across the interface between the two fluids for viscous fingering.10,31,32 For Saffman–Taylor viscous fingering, the flow in each of the fluid phases is governed by Darcy's law33,34
v = κ μ p ,
(3)
where v is the fluid velocity, while κ and μ are, respectively, the permeability of the porous medium and the fluid's viscosity. Taking into account the flow incompressibility, · v = 0 for incompressible fluids, the Laplace equation P 2 (P being the pressure field) holds everywhere in both fluids. It follows from this that the same Eq. (2) holds for the pressure drop across the interface, and thus, the interface displacement is driven by Laplacian growth if the fluids' viscosities are different.
The linearity between the pressure drop and the velocity or the flow rate indicated by Eq. (3) does not, however, hold for steady-state flow in a certain range of capillary numbers. Steady-state flow implies the simultaneous flow of both fluids in the porous medium for a sufficiently long time so that the system has reached a state when the statistical averages of macroscopic quantities, such as the saturation or the global pressure drop, do not drift with time anymore, while both fluids are still flowing. Experiments in a porous medium consisting of a monolayer of glass beads in which air and a water–glycerol mixture were displaced simultaneously have shown that the total flow rate Q in steady state varied with the applied pressure drop Δ P as a non-linear power law of exponent 1.85.35,36 Later, experiments with two incompressible fluids in the same porous medium found the exponent to be 1.35 or 1.5 depending on the fractional flow.37 Experiments with three-dimensional porous media made of glass bead packings38–40 and real core samples,41–43 performed by different groups, have further established this non-linearity with different exponent values. Various theoretical approaches35,36,44,45 and numerical modelings with variable-radii tubes,46 capillary bundles,47 pore-networks models,40,44 and lattice Boltzmann simulations48 have been carried out to understand the origin of the non-linearity that makes the rheology to deviate from Eq. (3). In general, the non-linear power law can be expressed by the expression
Q ( Δ P P t ) β ,
(4)
where P t is a threshold pressure that may exist so that there will be no flow in the system below P t, and β > 1 is the non-linear exponent. Fundamentally, this non-linearity is related to the distribution of the capillary barriers at the pore throats, which result from surface tension at the menisci between the two fluids; disorder in the pore geometry induces disorder in the spatial distribution of the capillary barriers. With the increase in Δ P, more and more pores progressively allow for the barrier to be overcome and flow to be conducted, which makes Q increase faster than the increase in Δ P, and hence, β > 1.49 When all the available pores along the interface, at any time, start flowing at a sufficiently high Δ P, the relation becomes linear. For a simplified porous system such as a capillary bundle, the value of β can be determined analytically by integrating the individual flow rates of the tubes over the distribution of the capillary barriers.47 For a porous network, the distribution of pore throat sizes50 and the wettabilities51,52 strongly control the value of β, as the distribution of capillary barriers depends on them.

In this article, we investigate to which extent this non-linearity in the flow rate, which has been demonstrated to exist for the steady state due to disorder in the capillary barriers, also exists for the growth of viscously-unstable fingers during drainage, and what its characteristics are. Equation (2) for DLA considers p = 0 at the perimeter of the aggregate, which, fingering in two-phase flow, would correspond to the absence of surface tension. The analogy to DLA therefore only applies to capillary numbers that are sufficiently high for capillary forces to be insignificant as compared to viscous forces. In the intermediate regime between capillary and viscous fingerings, we may expect that capillary forces will compete with viscous forces, so that the disorder in capillary barriers will play a role in the finger growth. Løvoll et al.10 and Toussaint et al.53 have investigated this displacement regime experimentally using two-dimensional (2D) porous media consisting of Hele–Shaw cells filled with a monolayer of glass beads, and measured the fractal dimension of the fingers and the statistical width of the front,29,54,55 comparable to λ in Eq. (1). Based on these measurements, they showed that the drainage fingers are characteristically different from DLA. They further proposed a quadratic relationship between the finger growth rate and the local pressure gradient in the intermediate flow regime, resulting from the pore-scale disorder. For DLA however, an ad hoc-type surface tension was considered by introducing a sticking probability,31 which showed no change in the fractal dimension.56,57 It was then argued that the two-phase flow fingers are more similar to dielectric breakdown models (DBM) (or the η models),58–60 which are generalizations of DLAs with power-law relationships of exponent η between the growth probability and the local growth/displacement-driving field.

The experimental study on 2D porous media10,53 left a number of open questions, in particular concerning the dependence of the non-linear growth exponent and statistical width of the viscous fingers on the capillary number. In this paper, we study viscously-unstable drainage as a function of the capillary number using large-scale simulations in dynamic pore networks. We inject a lower-viscosity fluid at one edge of the network filled with a higher viscosity fluid and characterize the statistical properties of the resulting viscous finger in a reference frame attached to the most advanced fingertip, as the finger grows inside the porous medium. We explore in particular the relationship between the finger's growth rate and the local pressure drop across the interface, and how it depends on the capillary number. We also focus on the statistical width parameter λ, analog to the λ in Eq. (1), compare it with the values obtained in the previous experimental study as well as in studies of DLA and Saffman–Taylor viscous fingering, and characterize its dependence on the capillary number.

In the following, we first describe the computational model and the simulation procedure in Sec. II, which is based on dynamic pore-network modeling61 with a specific boundary condition. In Sec. III, we present our results, where we first characterize the statistical profiles related to the volume and growth of the fingers in Subsection III A. In Subsection III B, we characterize the shape profile and measure the width ratio λ. In Subsection III C, we then explore the relationship between the growth rate and the local pressure drop across the interface. Finally, we discuss the results and provide an overall conclusion in Sec. IV.

The core of our simulation consists of a dynamic pore-network model,61 which we tweaked to adapt to the present problem. The model has been developed for over a decade62 and tested against numerous experimental,40,63 theoretical,44,64–66 and lattice Boltzmann simulations.67 In this computational method, the pore space of a porous medium is modeled by a network of links and nodes. Our network is spread in two dimensions (2D) within a flat cuboid domain similar to a Hele–Shaw cell and consists of N W × N L number of links embedded in a diamond lattice. Such a network with NW = 20 and NL = 80 is shown in Fig. 1. Here, the subscripts W and L refer to the directions of the network orthogonal and parallel to the overall flow, respectively. Each link of the network is of length l = 1 mm. The total network is therefore W = h N W mm wide and L = h N L mm long, where h = l / 2. We assign the entire pore space of the network to the links, so the nodes only represent the position of the link intersections. The correspondence between this geometry and that of a granular quasi-two-dimensional porous medium consisting of cylinders is illustrated in Fig. 2(a), where a small part of the porous medium is shown. The links are therefore composite in nature, which means that each link must represent a narrow pore throat in between two wider pore bodies. Such a geometry is modeled by an hourglass-shaped converging–diverging tube [as illustrated in Fig. 2(b)] whose cross-sectional area varies along its length. This results in a variation in the capillary pressure p as the interface (i.e., meniscus) between the two fluids moves along a link. We model this variation with a modified Young–Laplace equation45,62,68
| p i ( ϵ ) | = 2 γ r i [ 1 cos ( 2 π ϵ l ) ] ,
(5)
where y is the position of an interface inside a link i. Such variation in the capillary pressure is shown in Fig. 2(b). Here, γ = γ cos θ, where γ and θ represent the surface tension between the two fluids and the contact angle, respectively. The ri is the average radius of the ith link, the value of which we choose from a uniform distribution of random numbers in the range between 0.1 l and 0.4 l, which introduces the disorder in the network.
FIG. 1.

A pore-network of dimension L × W consisting of NW = 20 and NL = 80 links. The length of each link is l, and therefore, h = l / 2. The network was initially filled with high-viscosity fluid (represented in gray). The low-viscosity fluid (represented in blue) is injected with a total flow rate Q through NI = 2 inlet nodes at the bottom, indicated by the black dots. All the other nodes at the bottom boundary and all on the two vertical boundaries are closed, which is indicated by the thick black lines. The nodes at the top boundary are open and work as outlets. The (x, z) coordinate axes indicated by the orange lines are attached to the most advanced fingertip and move with it. The part of the network between the two red horizontal lines is used for the measurements of different quantities in order to avoid finite-size effects at the inlet and the outlet.

FIG. 1.

A pore-network of dimension L × W consisting of NW = 20 and NL = 80 links. The length of each link is l, and therefore, h = l / 2. The network was initially filled with high-viscosity fluid (represented in gray). The low-viscosity fluid (represented in blue) is injected with a total flow rate Q through NI = 2 inlet nodes at the bottom, indicated by the black dots. All the other nodes at the bottom boundary and all on the two vertical boundaries are closed, which is indicated by the thick black lines. The nodes at the top boundary are open and work as outlets. The (x, z) coordinate axes indicated by the orange lines are attached to the most advanced fingertip and move with it. The part of the network between the two red horizontal lines is used for the measurements of different quantities in order to avoid finite-size effects at the inlet and the outlet.

Close modal
FIG. 2.

(a) Illustration of the network of pores and nodes. The pore space is indicated by white color, whereas the gray circles represent the grains in a porous rock or the glass beads in a Hele–Shaw cell. One of the links is colored by dark gray. The links intersect at nodes at the intersections of the dashed lines. (b) Shape of individual pores in the mid-horizontal plane of the network, and variation of the capillary pressure p , given by Eq. (5), along the length of the pore. The white and gray segments inside the pore represent, respectively, the wetting and non-wetting fluid blobs; there are three interfaces between them in this case, which are in contact with the beads at positions ϵ1, ϵ2, and ϵ3 along the pore length.

FIG. 2.

(a) Illustration of the network of pores and nodes. The pore space is indicated by white color, whereas the gray circles represent the grains in a porous rock or the glass beads in a Hele–Shaw cell. One of the links is colored by dark gray. The links intersect at nodes at the intersections of the dashed lines. (b) Shape of individual pores in the mid-horizontal plane of the network, and variation of the capillary pressure p , given by Eq. (5), along the length of the pore. The white and gray segments inside the pore represent, respectively, the wetting and non-wetting fluid blobs; there are three interfaces between them in this case, which are in contact with the beads at positions ϵ1, ϵ2, and ϵ3 along the pore length.

Close modal
The capillary numbers for both the viscous and capillary fingering regimes are sufficiently small for any inertial effect to be negligible. For fully developed laminar flow of incompressible fluids in each link, we therefore consider the following equation for the flow rate of the fluids:23,69
q i = k i l μ i [ Δ p i b p i ( ϵ b ) ] ,
(6)
where Δ p i is the pressure drop between the two nodes across the ith link. The network is assumed to be placed horizontally, and no gravity is considered. Here, μi is the effective viscosity of the two fluids inside the link, which is given by μ i = s i μ n + ( 1 s i ) μ w, where μ n and μ w are the viscosities of the non-wetting and wetting fluids, respectively, and si is the non-wetting saturation in the ith link, i.e., the proportion of the link volume that is occupied by the non-wetting fluid. The term ki represents the mobility of the link, given by k i = a i r i 2 / 8 for the Hagen–Poiseuille flow in a circular cross section, where a i = π r i 2 is the cross-sectional area of the link. The summation of the capillary pressures p i in Eq. (6) is over all the interfaces inside the link i, obtained using Eq. (5). For the link shown in Fig. 2(b), for example, the summation will be over the three interfaces at ϵ1, ϵ2, and ϵ3.

To find out the pressures pj at every node j of the network, we use the Kirchhoff equation for incompressible flow i n j q i = 0 for each node. Here the summation is over the links nj connected to a node j. This provides a set of linear equations which we solve by conjugate-gradient method70 with proper boundary conditions. This is illustrated in Fig. 1 where the two nodes (NI = 2) at the bottom edge act as the inlets. All the nodes at the top edge work as the outlets through which fluids leave the system. The two vertical boundaries of the network parallel to the overall flow are closed. The network is initially filled with wetting fluid (gray colored in the figures), and we inject the non-wetting fluid (blue colored in the figures) through inlet nodes with a total constant flow rate Q. We therefore set pj = 0 for all the outlet nodes and q i = Q / N I for all the virtual links connected to the inlet nodes, indicated by the two arrows in Fig. 1. These serve as the boundary conditions for solving the linear system of equations for the node pressures.

The positions of the two fluids inside every link are assigned by the positions of the interfaces between the blobs of the two immiscible fluids in every link as shown in Fig. 2(b). At every time step, the displacements of the two fluids are performed by updating these positions by a distance
Δ ϵ i = Δ t q i / a i
(7)
in the direction of the flow in the corresponding link. Here, Δ t is a time step chosen in such a way that the displacement of an interface inside any link does not exceed more than 0.1 l.

One final detail about the modeling is how to distribute the fluids from links to their neighboring links through the nodes. For this, first, the links carrying fluids toward and away from every node are identified, we call them the incoming and outgoing links, respectively. Then, for every node j, the total volumes of the wetting and non-wetting fluids ( V j w and V j n, respectively, and V j = V j w + V j n) arriving from the incoming links to the node are measured using Eq. (7). These volumes are then distributed to each outgoing links i by placing new wetting and non-wetting blobs of volumes V i w = q i Δ t V j w / V j and V i n = q i Δ t V j n / V j, respectively. This means that the ratio of V i w to V i n in any outgoing link is identical to the ratio between the incoming wetting and non-wetting volumes at the corresponding node, and the ratio between the total volumes ( V i w + V i n) injected in different outgoing links i is the same as the ratio between the flow rates qi in those links. Furthermore, technical details about this algorithm can be found in Ref. 61.

We performed simulations of drainage displacement at constant flow rate Q set by the capillary number Ca, which is a dimensionless number that quantifies the typical ratio of the viscous to the capillary forces across a pore positioned at the fluid–fluid interface. For a single tube and if one of the fluids' viscosity is negligible in comparison with the other fluid's viscosity, Ca is expressed as12,71
Ca = μ v γ ,
(8)
where v is the velocity of the fluids in the tube and μ is the viscosity of the high-viscosity fluid. For a porous medium consisting of many pores and without any assumption on the viscosity, this expression was extended by many authors where the viscosity μ was referred to as the viscosity of the injected fluid into the porous medium.23,72 However, in the case of a low-viscosity fluid injected into a fluid of much higher viscosity, as is the case in the configurations which we consider, the viscous forces are controlled by the higher viscosity, and μ in Eq. (8) is therefore considered to be the viscosity of the defending fluid in such cases.10,73 In experiments, the velocity v in Eq. (8) for a porous medium is usually considered as the Darcy or filtration velocity of the invading phases, which is the total flow rate per unit cross-sectional area of the porous medium. However, for our pore network, only the areas of the channels are properly defined, and therefore, we define v as the true mean longitudinal fluid velocity, equal to the flow rate per unit pore area. Furthermore, Eq. (8) assumes that the permeability can be equated to the square of the typical channel cross section, whereas many experimental studies including those addressing the same topic,10,53 also considered that the viscous pressure drop should involve the medium's permeability, which leads to further differences in the definition of Ca. In addition, the wetting angle could also be taken into account in the capillary pressure drop, thus appearing in the definition of Ca as well. When comparing our simulation results to experimental results under identical flow conditions, the similarity of the flow conditions must be assessed based on Ca estimates obtained from the same definition, so a conversion of the Ca values (either experimental or numerical) may be necessary.
In order to find out the value of Q that must be imposed in the simulation for a given Ca, we first define it as74 
Ca = Δ p visc Δ p cap = Δ P W h / L 2 γ / r ¯ ,
(9)
where Δ p visc = Δ P W h / L and Δ p cap = 2 γ / r ¯ are, respectively, the average viscous and capillary pressure drops across a pore. Here, Δ P W is the total viscous pressure drop across the network when it is completely filled with the wetting fluid and r ¯ is the average link radius. We then first solve the viscous pressure drop Δ P visc , 1 for a test flow rate Q1 for the network saturated with wetting fluid, and then determine the flow rate Q by using
Q = 2 γ L h r ¯ Q 1 Δ P visc , 1 Ca .
(10)
This is because Q Δ P W for single-phase flow. The simulation is then performed with q i = Q / N I for all the virtual inlet links.

We have considered a network of dimension 32 × 160 links in this study, and this is based on the maximum computational time we could spend reasonably. We had to make additional sweeps through the entire network at every time step to compute the volume density, the growth density, and the local pressure drops, which we will discuss later. This made the simulations considerably more time-consuming compared to usual pore-network simulations with the same model. Furthermore, in the analyses, we discarded the regions of the network that are distant from the inlet and the outlet by less than NW rows, as shown by the red lines in Fig. 1, in order to avoid any boundary effect. In Table I, we show the numbers of different realizations of the network considered for different values of M and Ca. All the measured quantities in different plots in this article are averaged over these many samples. The surface tension between the fluids was chosen to γ = 0.03 N / m. Two viscosity ratios, M = μ n / μ w = 10 4 with μ n = 10 5 Pa s and μ w = 10 1 Pa s, and M = 10 5 with μ n = 10 5 Pa s and μ w = 1 Pa s, were considered. A typical set of simulations for different capillary numbers and viscosity ratios are shown in Fig. 3, where the black and blue colors represent the wetting and non-wetting fluids, respectively. The direction of the overall flow in these images is from the bottom to the top as indicated by the arrow. The two red lines indicate the aforementioned regions near the inlets and the outlets that are not taken into account in the analyses. The gray shades indicate the average volume density, which we describe in the following.

TABLE I.

Number of samples simulated for different values of Ca and M. The values of M are 10 4 and 10 5, whereas the values of Ca are in the range 0.01–0.9. All the measured quantities are statistically averaged over these many samples for the respective values of Ca and M.

Ca 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90
M = 10 4  57  79  93  105  112  117  123  125  125  125  121  105  98  90  84  81  76  73 
M = 10 5  30  49  56  66  71  74  81  81  82  81  80  65  62  57  54  51  46  45 
Ca 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90
M = 10 4  57  79  93  105  112  117  123  125  125  125  121  105  98  90  84  81  76  73 
M = 10 5  30  49  56  66  71  74  81  81  82  81  80  65  62  57  54  51  46  45 
FIG. 3.

Drainage fingers generated for viscosity ratio M = 10 4, and for different values of Ca. The blue and black colors, respectively, denote the invading non-wetting and the defending wetting fluids. The inlets are at the bottom edge of the networks, and the overall flow direction is toward the top. The gray shade indicates the non-wetting volume density σ ̃ ( x , z ), averaged over time and samples in the moving coordinate system.

FIG. 3.

Drainage fingers generated for viscosity ratio M = 10 4, and for different values of Ca. The blue and black colors, respectively, denote the invading non-wetting and the defending wetting fluids. The inlets are at the bottom edge of the networks, and the overall flow direction is toward the top. The gray shade indicates the non-wetting volume density σ ̃ ( x , z ), averaged over time and samples in the moving coordinate system.

Close modal
The volume density of the fingers represents a statistical map of the occupation of the two-dimensional network space by the invading non-wetting fluid. This provides a statistical shape of the displacement structure, which is comparable with the smooth continuum-medium fingers.12 The density maps are functions of the distance from the most advanced fingertip, z, as the displacement process can be considered statistically stationary in an (x, z) reference frame attached to that tip (see Fig. 1), as previously shown for experimental drainage fingers,10,53 as well as for DLA fingers.29,54 With respect to this dynamic (x, z) reference frame, we define the volume density σ ̃ ( x , z ) as
σ ̃ ( x , z ) = 1 h 2 n ( x , z ) k = 1 n ( x , z ) V k ( x , z ) ,
(11)
which is essentially the volume of injected fluid per unit area at a given position in the (x, z) reference plane, averaged over time, V k ( x , z ) being the non-wetting volume in the link positioned at position (x, z) at time step k. Here, we point out that the lowest unit of resolution for our measurements is one single link, which means that we have only one data point within an area of h2. Furthermore, as the available window of the network for the measurements varies with time, the number of elements within the summation is different for different links. This is accounted by n(x, z) in the above expression, which is the number of times the volume in the link at (x, z) is added in the summation. This also means that the data near the fingertip contain better statistical averages compared to those far behind the tip. The map is then further averaged over a number of different samples of the network as indicated in Table I to obtain statistically averaged σ ̃ ( x , z ) map. In Fig. 3, we show the spatial distribution of σ ̃ ( x , z ) as gray-scale maps where the dark and light shades correspond to the smallest and largest values, respectively. For an observer sitting at the origin of the moving coordinate system, this volume profile of the finger will remain the same statistically during the invasion of the fingers in an infinitely long network, and only the length of the profile will increase in the z direction.
We next measure the volumetric growth rate ν ̃ ( x , z ) of the invading fingers, which is defined as the average increase in the volume of the finger per unit time and unit area at (x, z),
ν ̃ ( x , z ) = 1 h 2 n ( x , z ) k = 1 n ( x , z ) 1 Δ t k [ V k ( x , z ) V k 1 ( x , z ) ] ,
(12)
where Δ t k is the time interval between the time steps k and k – 1. Note that ν ̃ ( x , z ) has the dimension of a velocity. By integrating ν ̃ ( x , z ) in both x and z directions, we have
0 L 0 W ν ̃ ( x , z ) d x d z = Q
(13)
as the total increase in the volume of the whole finger per unit time is equal to the volumetric rate of the steady fluid injection into the system. By normalizing ν ̃ ( x , z ) by Q, we then define the growth density as
ϕ ̃ ( x , z ) = 1 Q ν ̃ ( x , z ) .
(14)
We assume that the system is statistically symmetric with respect to the longitudinal mid-section ( x = w / 2) of the flow cell, and we define the quantities σ ( z ) , ν ( z ), and ϕ ( z ) in the longitudinal direction as a function of the sole distance z from the tip, by integrating σ ̃ ( x , z ) , ν ̃ ( x , z ), and ϕ ̃ ( x , z ) in the x direction,
σ ( z ) = 0 W σ ̃ ( x , z ) d x , ν ( z ) = 0 W ν ̃ ( x , z ) d x , ϕ ( z ) = 0 W ϕ ̃ ( x , z ) d x
(15)
from which we can also define ϕ ( z ) = ν ( z ) / Q. In Fig. 4, we plot h ϕ ( z ) as a function of the number of rows from the most advanced tip, z/h, for different values of Ca and M. Near the fingertips, the data show an exponential-type decay of ϕ ( z ). This is indicated by the solid line where we plot a model function B exp ( z h ξ ) with ξ = 4.6 and B = 0.18. Here, ξ is a characteristic decay length or screening length, which characterizes the active invasion zone. The value of ξ appears to be the same for the whole range of capillary numbers and viscosity ratios, showing that it is a characteristic of the porous medium, which was also observed in the drainage experiments Hele–Shaw cell filled with a monolayer of glass beads.10 The same data are plotted in a linear-logarithmic scale in the insets, which show that the exponential decay is valid for z < W / 2 for lower capillary numbers, whereas the validity increases up to z < W at high capillary numbers. This is indicated by the linear dashed line, which has a slope 1 / ξ and intercept ln ( B ). This characteristic is also in agreement with experiments,10 which showed similar exponential decay of finger growth within the same range from the tip. Away from the tip, ϕ ( z ) is negligibly small as seen from the scattered data points in the insets, which means that the active growth of the fingers happens near the fingertip, whereas they are almost frozen far behind the tip.
FIG. 4.

Plot of the dimensionless growth density h ϕ ( z ) as a function of the row numbers from the tip, z/h, for the viscosity ratios (a) M = 10 4 and (b) M = 10 5. The different sets in each plot correspond to different values of Ca indicated in the legend of both plots. The solid line corresponds to the model function B exp ( z h ξ ) with ξ = 4.6 and B = 0.18, for both plots. The inset shows a linear-logarithmic plot of the same data where the dashed line has a slope of 1 / ξ and intercept of ln B.

FIG. 4.

Plot of the dimensionless growth density h ϕ ( z ) as a function of the row numbers from the tip, z/h, for the viscosity ratios (a) M = 10 4 and (b) M = 10 5. The different sets in each plot correspond to different values of Ca indicated in the legend of both plots. The solid line corresponds to the model function B exp ( z h ξ ) with ξ = 4.6 and B = 0.18, for both plots. The inset shows a linear-logarithmic plot of the same data where the dashed line has a slope of 1 / ξ and intercept of ln B.

Close modal
Note that σ ( z ) and ϕ ( z ) are related to each other. Indeed, let us now examine the displacement process in the referential ( x , z ̃ ) attached to the inlet boundary of the network, with z ̃ and z related to each other according to z = z ̃ tip ( t ) z ̃ = v tip t z ̃, where v tip is the velocity of the most advanced tip of the non-wetting fluid phase and z ̃ tip ( t ) is its position along z ̃. That velocity is constant in time, as shown in Fig. 5, where z ̃ tip is plotted as a function of the time normalize by the time t scaled by tf, where tf is the time at which the tip reaches the outlet boundary of the porous network. At any time t, the volume of non-wetting fluid in a slice of the network normal to the mean flow direction and located between z ̃ and z ̃ + Δ z ̃ is σ ̂ ( z ̃ , t ) Δ z ̃, with Δ z ̃ = Δ z. It results from the cumulative growth of the fingers between the time t 0 = z ̃ / v tip at which the most advanced tip first arrived at z ̃ and t, and can thus be expressed as
σ ̂ ( z ̃ , t ) Δ z = Δ z z ̃ / v tip t ν ̂ ( z ̃ , u ) d u .
(16)
The steady volumetric density defined above in the (x, z) reference frame, σ ( z ), is related to σ ̂ ( z ̃ , t ) through the following relation:
t σ ( z ) = σ ̂ ( v tip t z , t ) .
(17)
A similar expression relates the volumetric growth rate expressed in the moving reference frame, ν ( z ), and that expressed in the laboratory frame, ν ̂ ( z ̃ , t ). Hence,
t σ ( z ) Δ z = Δ z t 0 t ν ̂ ( v tip t z , t ) d t , i . e . , σ ( z ) Δ z = Δ z 0 z ν ( z ) d z v tip , i . e . , σ ( z ) Δ z = Q Δ z v tip 0 z ϕ ( z ) d z .
(18)
Defining the cumulative growth density in the z direction as Φ ( z ) = 0 z ϕ ( z ) d z and σ c = Q / v tip leads to
Φ ( z ) = σ ( z ) σ c ,
(19)
which relates the volume density with the cumulative growth density.
FIG. 5.

Plot of the position of the most advanced fingertip, z ̃ tip, as a function of the scaled time t / t f, where t f is the time when the most advanced fingertip reaches the final position in the simulation. The two plots correspond to the viscosity ratios (a) M = 10 4 and (b) M = 10 5; the different symbols correspond to the different values of Ca indicated in the plot legends.

FIG. 5.

Plot of the position of the most advanced fingertip, z ̃ tip, as a function of the scaled time t / t f, where t f is the time when the most advanced fingertip reaches the final position in the simulation. The two plots correspond to the viscosity ratios (a) M = 10 4 and (b) M = 10 5; the different symbols correspond to the different values of Ca indicated in the plot legends.

Close modal

To verify this functional form, the tip velocity is first measured from the slopes as the straight lines in Fig. 5 multiplied by tf. We then plot σ ( z ) scaled by σc, for different values of Ca and M (Fig. 6). Notice that, though σ ( z ) is a function of the capillary numbers as shown in the insets, all σ ( z ) / σ c plots collapse on a single master curve for different values of Ca. With Eq. (19), it therefore implies that there is one unique cumulative growth density Φ ( z ), independent of Ca. We therefore calculate the average cumulative growth density Φ ( z ) from the data in Fig. 4 and plot it in Fig. 6 as a solid black line. The line shows an excellent agreement with the σ ( z ) / σ c plots, as expected from Eq. (19). Such a growth property of the front is characteristic of the viscous-dominated regime, or of a regime where viscous forces compete with capillary forces, and therefore will not be observed in the pure capillary fingering regime where the invasion is completely controlled by the disorder in the capillary thresholds.

FIG. 6.

Plot of the scaled longitudinal volume density σ ( z ) / σ c as a function of the row numbers z/h behind the most advanced fingertip. The two figures correspond to the viscosity ratios (a) M = 10 4 and (b) M = 10 5, and the different data sets correspond to simulations with different injection rates Q, and thus to different capillary numbers Ca. The values of Ca are indicated in the figure legends. The average cumulative growth density Φ ( z ) is plotted as solid black line, which shows the agreement with Eq. (19). The black, blue, and red dashed lines are plotted using Eq. (1) with λ = 0.4, 0.5, and 0.62, which correspond to the experimental observation for porous Hele–Shaw cell,10,53 the DLA,55,75 and the Saffman–Taylor solution for low surface tension limit,12 respectively. The insets show the variations of unscaled σ ( z ) for different values of Ca.

FIG. 6.

Plot of the scaled longitudinal volume density σ ( z ) / σ c as a function of the row numbers z/h behind the most advanced fingertip. The two figures correspond to the viscosity ratios (a) M = 10 4 and (b) M = 10 5, and the different data sets correspond to simulations with different injection rates Q, and thus to different capillary numbers Ca. The values of Ca are indicated in the figure legends. The average cumulative growth density Φ ( z ) is plotted as solid black line, which shows the agreement with Eq. (19). The black, blue, and red dashed lines are plotted using Eq. (1) with λ = 0.4, 0.5, and 0.62, which correspond to the experimental observation for porous Hele–Shaw cell,10,53 the DLA,55,75 and the Saffman–Taylor solution for low surface tension limit,12 respectively. The insets show the variations of unscaled σ ( z ) for different values of Ca.

Close modal

The average volume densities of the fractal fingers can be used to compare their shapes statistically with the smooth shapes of continuum-medium fingers, as well as with the shapes of DLA. By mapping the volume density function σ ( z ) to the Saffman–Taylor equation [Eq. (1)] for continuum viscous fingers, one can calculate the width ratio λ and use that to compare the statistical shape profile of different fingers. Before presenting the detailed measurements of λ, we show a quick comparison in Fig. 7 where only the part of the network with volume density σ ̃ ( x , z ) σ max / 2 are shown. This region corresponds to the porous media with large invasion54,55,76 and is comparable to the Saffman–Taylor equation. For the viscous fingers in the continuum Hele–Shaw cell without any porous structure, the width ratio λ in Eq. (1) has a value of 0.5 in the limit of negligible surface tension, and it increases with the increase in surface tension.21,77 For experiments with porous Hele–Shaw cells on the other hand, λ was found to be around 0.4,10,53 whereas for off-lattice DLAs, λ was found to be around 0.62 for linear channels.55,75 The solutions of Eq. (1) for these three values are shown in Fig. 7, where the green, yellow, and blue lines correspond to λ = 0.4, 0.5, and 0.62, respectively. The lines show deviations from the gray-scale region, and the deviations also depend on Ca. We also show the solutions of Eq. (1) for these three values of λ in Fig. 6 by the dashed lines, which again show some deviations from the results for the drainage fingers here.

FIG. 7.

Volume density maps thresholded with σ ̃ ( x , z ) σ ̃ max / 2 for M = 10 4 and for Ca = 0.04, 0.08, and 0.1, as indicated below the figures. The gray shades correspond to the same density map as shown in Fig. 3 for the respective values of Ca, but now only the links preserved by the cutoff are shown. The green, yellow, and blue lines are drawn using the Saffman–Taylor equation [Eq. (1)] with λ = 0.4, 0.5 and 0.62, which correspond to the experimental observation for porous Hele–Shaw cell, the continuum-medium Saffman–Taylor finger and the DLA finger, respectively.

FIG. 7.

Volume density maps thresholded with σ ̃ ( x , z ) σ ̃ max / 2 for M = 10 4 and for Ca = 0.04, 0.08, and 0.1, as indicated below the figures. The gray shades correspond to the same density map as shown in Fig. 3 for the respective values of Ca, but now only the links preserved by the cutoff are shown. The green, yellow, and blue lines are drawn using the Saffman–Taylor equation [Eq. (1)] with λ = 0.4, 0.5 and 0.62, which correspond to the experimental observation for porous Hele–Shaw cell, the continuum-medium Saffman–Taylor finger and the DLA finger, respectively.

Close modal
In order to measure the value of the width parameter λ, we measure a transverse volume density ρ ( x ) in a region where the growth is almost frozen, that is, far behind the most advanced fingertip. If L f is the length of this frozen zone in the z direction, then one has
0 W W W + L f σ ̃ ( x , y ) d z d x = σ c L f
(20)
as L f increases with the same velocity v tip. We therefore define ρ ( x ) as
ρ ( x ) = W σ c L f W W + L f σ ̃ ( x , z ) d z
(21)
with 0 W ρ ( x ) d x = W. In Figs. 8(a) and 8(b), we plot ρ ( x ) as a function of the scaled transverse position x/W for the two different viscosity ratios. The datasets show a maximum at the middle and then decrease on both sides. This is an expected behavior, similar to what is observed for DLA55 and for the experiments of viscous fingers in porous media.10,53 However, the plots drift with Ca, indicating a variation of λ with Ca. From these data, the width ratio λ can be estimated in two ways: either from (A) λ = ( x + x ) / W, where ρ ( x + ) = ρ ( x ) = λ max / 2 or alternatively from (B) λ = 1 / ρ max.54,76 The two measurements are shown in Fig. 8(c), where λ seems to vary systematically with Ca. Interestingly, the lowest value here is at 0.5, which is the solution for the continuum Saffman–Taylor fingers for negligible surface tension limit.12 Higher values of λ were found for increasing surface tension for the continuum fingers.21,77,78 The value of λ 0.62 for DLAs in linear channels55,75 coincides with some of the data points here; however, the two-phase flow experiments in porous Hele–Shaw cells showed a much lower value of λ 0.4.10,53 Note however that in contrast to the present system, the capillary threshold distribution was not uniform in these experiments, despite the fact that a uniform distribution was assumed for some of the interpretations. Here, for our system, λ increases up to a much higher value of around 0.8. Such values are similar to the fingers in divergent cells, which are linear cells with an angular wedge at the inlet, a geometry that is in between the linear and circular geometries.79 There, λ was found to increase when increasing the angle of the wedge, for example, 0.77–0.82 for a wedge angle of 90°.29,55 Note that, in our simulations, we inject fluid only through a few of the nodes at the center of the network's inlet boundary, and all other nodes on the two sides of these injection nodes are blocked, see Fig. 1. This may have a wedge effect on the value of λ, similar to the divergent inlets. However, the fact that here λ shows a minimum at an intermediate Ca and increases on both sides indicates that there is a combined effect of surface tension and inlet geometry, and it therefore needs further in-depth study. We leave this for the future, as our main focus here is to explore the effect of disorder on the relationship between the growth and pressure drop.
FIG. 8.

Plots of the average transverse volume density profile ρ ( x ) as a function of the transverse normalized coordinate x/W are shown for (a) M = 10 4 and (b) M = 10 5. The different symbols show the results for different values of Ca, and the solid line shows the profile averaged over all datasets. The values of λ measured using the two different methods, (A) λ = 1 / ρ max and (B) λ = ( x + x ) / W, where ρ ( x + ) = ρ ( x ) = λ max / 2, are shown in (c) as a function of the capillary number.

FIG. 8.

Plots of the average transverse volume density profile ρ ( x ) as a function of the transverse normalized coordinate x/W are shown for (a) M = 10 4 and (b) M = 10 5. The different symbols show the results for different values of Ca, and the solid line shows the profile averaged over all datasets. The values of λ measured using the two different methods, (A) λ = 1 / ρ max and (B) λ = ( x + x ) / W, where ρ ( x + ) = ρ ( x ) = λ max / 2, are shown in (c) as a function of the capillary number.

Close modal
The volumetric growth of the region occupied by the invading fluid inside a pore depends on whether the viscous pressure drop between the invading and defending fluids is able to overcome the capillary barrier inside the pore and thus to displace the interface between the two fluids. Assuming that the movement of the interfaces follows Eq. (6), a theoretical methodology was suggested in Refs. 10 and 53 to find the relationship between the growth rate and the local pressure drop across the fluid–fluid interface by integrating over the distribution of the capillary barriers. We follow their approach in the following and modify it to adapt our system. If P is the threshold pressure required to invade a link and g ( P ) is the normalized distribution of the thresholds over the network, we may integrate over g ( P ) to obtain the growth rate
ν ̃ ( x , z ) = 1 h 2 q ( x , z ) g ( P ) d P = κ h 2 ( Δ p ( x , z ) P ) Θ ( Δ p ( x , z ) P ) g ( P ) d P ,
(22)
where Δ p ( x , z ) is the pressure drop between the invading and defending fluids inside a link i at (x, z) and q(x, z) is the corresponding flow rate of the fluids in that link. Here, we consider κ = k ¯ / μ w, where k ¯ is the average mobility of the pores and μw is the viscosity of the high-viscosity defending fluid. A pore will be invaded only when Δ p ( x , z ) > P , hence the Heaviside function Θ ( Δ p ( x , z ) P ).
If we consider P to be uniformly distributed between P t and P m, in the form
g ( P ) = 1 G Θ ( P P t ) Θ ( P m P ) ,
(23)
where G = P m P t, then Eq. (22) becomes
ν ̃ ( x , z ) = κ h 2 G P t P m ( Δ p P ) Θ ( Δ p P ) d P ,
(24)
where Δ p Δ p ( x , z ). If the viscous pressure drop exceeds the capillary threshold for all pores, i.e., if Δ p ( x , z ) > P m, then Eq. (24) reduces to
ν ̃ ( x , z ) = κ h 2 G P t P m ( Δ p P ) d P = κ 2 h 2 ( P m P t ) [ ( Δ p P ) 2 ] P t P m , ν ̃ ( x , z ) = κ h 2 [ Δ p P m + P t 2 ] .
(25)
The growth rate in this regime then varies linearly with the excess pressure drop with respect to the mean capillary threshold ( P m + P t ) / 2.
For Δ p ( x , z ) < P t on the other hand, Eq. (24) reduces to
ν ̃ ( x , z ) = κ h 2 G P t Δ p ( Δ p P ) d P , ν ̃ ( x , z ) = κ 2 h 2 ( Δ p P t ) 2 P m P t ( Δ p P t ) 2 .
(26)
The growth rate thus varies quadratically with [ Δ p ( x , z ) P t ]. This regime corresponds to moderate capillary numbers where the pressure drop Δ p ( x , z ) competes with the capillary thresholds, and the number of invaded pores thus increases with an increase in Δ p ( x , z ).
If we further assume that Δ p ( x , z ) is a function of z only, we may average Δ p ( x , z ) over x, which yields
ν ( z ) ( Δ P ( z ) P t ) 2 ,
(27)
where Δ P ( z ) = Δ p ( x , z ) x. Hence, in this regime, the average growth rate ν ( z ) depends quadratically on the excess pressure drop at z, [ Δ P ( z ) P t ]. Such a quadratic relationship between excess pressure drop and flow rate was also suggested for steady-state two-phase flow in porous media by mean-field calculations44 and observed by numerical simulations.40,44,48
Note however that the threshold distribution function g ( P t ) considered in the above derivation was assumed to be uniform. In our pore network, it is a more complex function, which is not straightforward to determine. It not only depends on the pore properties such as the size distribution50 and the pore wettabilities51,52 but also on the structure of the finger itself, that is, the distribution of the number of fingertips at any z. We therefore may assume a generalization of Eq. (27),
ν ( z ) ( Δ P ( z ) P t ) β ,
(28)
where β is an exponent of non-linearity. As mentioned in the Introduction, such non-linear relationships with β > 1 have been widely observed for steady-state flow in a certain regime of capillary numbers.35–48 
In order to examine the form of Eq. (28), we need to measure two quantities: (a) the average pressure difference between the defending and invading fluids, Δ P ( z ), as a function of z, and (b) the threshold pressure P t. At every time step, we obtain the pressures at the nodes of the network from the conjugate-gradient solver by solving Eqs. (6) and (5). An example is shown in Fig. 9, where the pressure values are indicated by different colors. Notice that the nodes connected by the invading non-wetting low-viscosity fluid (blue links) are at much higher pressures than those that are connected with the defending wetting fluid. The conjugate-gradient solver does not distinguish between the wetting and non-wetting pressures; we therefore identify, at every time step, the links filled with invading fluid as those with a non-wetting saturation larger than a threshold value ( s i > 0.98), and the nodes that are connected to them are marked as the non-wetting nodes. The rest of the nodes are then marked as the wetting nodes. We then measure the average excess wetting and non-wetting pressures as a function of the distance z from the most advanced fingertip as
P w ( z ) = 1 n w ( z ) k = 1 n w ( z ) p k ( z ) P 0 , and P n ( z ) = 1 n n ( z ) k = 1 n n ( z ) p k ( z ) P 0 ,
(29)
where n w ( z ) and n n ( z ) are the number of wetting and non-wetting nodes at z, and P0 is the average wetting pressure at the most advanced fingertip, z = 0. The pressure drop Δ P ( z ) between the wetting and non-wetting fluids at z is then calculated as
Δ P ( z ) = P n ( z ) P w ( z ) ,
(30)
which is then averaged over different time steps during the propagation of the fingers, and over different samples.
FIG. 9.

Node pressures (pj) obtained from conjugate-gradient solver, normalized by the maximum pressure p max at that time step, shown by colors for a network of 32 × 96 links. The color scale shows the normalized pressure values, where the darkest red corresponds to pj = 0 and the lightest yellow corresponds to the highest pressure. Here, Ca = 0.05 and M = 10 5. The injected fluid is colored blue. The box at the top right shows a close view of a small slice of the network. Notice that the nodes belonging to the low-viscosity finger are at much higher pressures than those belonging to the defending fluid.

FIG. 9.

Node pressures (pj) obtained from conjugate-gradient solver, normalized by the maximum pressure p max at that time step, shown by colors for a network of 32 × 96 links. The color scale shows the normalized pressure values, where the darkest red corresponds to pj = 0 and the lightest yellow corresponds to the highest pressure. Here, Ca = 0.05 and M = 10 5. The injected fluid is colored blue. The box at the top right shows a close view of a small slice of the network. Notice that the nodes belonging to the low-viscosity finger are at much higher pressures than those belonging to the defending fluid.

Close modal

In Fig. 10, P n ( z ) , P w ( z ), and Δ P ( z ) are plotted as functions of the row distance z/h from the fingertip for different values of Ca and M. We can notice a few things here. First, the wetting pressure P w shows a sharp variation with z, whereas the non-wetting pressure P n does not show any noticeable variation. This is because of the low viscosity of the invading fluid compared to the defending fluid. Second, the pressure drop Δ P ( z ) starts from a high value at z = 0, then falls rapidly within a distance z < W / 2 and almost goes to a plateau. This behavior is similar to the growth density ϕ ( z ) in Fig. 4, which showed a rapid exponential fall within a similar distance, indicating that the growth ϕ ( z ) is directly related to Δ P ( z ). Furthermore, we note that the plateau regimes of Δ P far behind the tip have finite values instead of zero. This however is the stagnant zone, where the fingers do not show any visible growth and ϕ ( z ) is essentially zero as seen in Fig. 4. Therefore, in this regime, the capillary barriers compete with the viscous pressure drop and prevent the interfaces from moving. The plateau value thus corresponds to the threshold pressure P t in Eq. (28), below which there is no growth. We can therefore measure P t from this plateau regime, by taking the average of Δ P ( z ) over its farthest data points. However, if we look carefully at the plots, we notice that not all data sets at the plateau exhibit a perfectly constant behavior, some of them are still decreasing slowly and therefore would need a much larger system to saturate. We therefore estimate the threshold pressure P t in one more way where we minimize the least-square fit error for Eq. (28). We try a set of trial values for P t around the plateau of Δ P ( z ) and then choose the one that provides the lowest least-square error for the straight-line fit for log [ ν ( z ) ] vs log [ Δ P ( z ) P t ]. These estimated values of P t are shown in Fig. 11 with blue triangles, where we also compare them with the measured values of P t, which are calculated as the average of Δ P over the farthest 20% rows from the fingertip. Notice that the values of P t from these two different calculations are very close to each other for the whole range of the capillary numbers. Furthermore, P t does not seem to depend on the viscosity ratio M, which confirms that the threshold is a quantity controlled solely by the capillary barriers, which are functions of the interfacial tension between the two fluids.

FIG. 10.

Plots of the non-wetting pressure P n ( z ), the non-wetting pressures P w ( z ), and the pressure drop Δ P ( z ), obtained using Eqs. (29) and (30), as a function of the distance from the most advanced fingertip z/h, for different capillary numbers. The top row corresponds to M = 10 4 and the bottom row to M = 10 5.

FIG. 10.

Plots of the non-wetting pressure P n ( z ), the non-wetting pressures P w ( z ), and the pressure drop Δ P ( z ), obtained using Eqs. (29) and (30), as a function of the distance from the most advanced fingertip z/h, for different capillary numbers. The top row corresponds to M = 10 4 and the bottom row to M = 10 5.

Close modal
FIG. 11.

Dependence of the threshold pressures P t on the capillary numbers Ca for the two different viscosity ratios (M), (a) 10 4 and (b) 10 5. The measured values of P t represented by the red circles are obtained by taking the averages of Δ P ( z ) for z / h > 76, i.e., over the last 20% of the data points on the plateau of Δ P ( z ) in Fig. 10. The estimated values of P t represented by the blue triangles are obtained by minimizing the least-square fit error of Eq. (28).

FIG. 11.

Dependence of the threshold pressures P t on the capillary numbers Ca for the two different viscosity ratios (M), (a) 10 4 and (b) 10 5. The measured values of P t represented by the red circles are obtained by taking the averages of Δ P ( z ) for z / h > 76, i.e., over the last 20% of the data points on the plateau of Δ P ( z ) in Fig. 10. The estimated values of P t represented by the blue triangles are obtained by minimizing the least-square fit error of Eq. (28).

Close modal

Finally, we now set out to verify the agreement with Eq. (28) and determine the exponent β. In Fig. 12, we plot the growth rate ν ( z ) as a function of the excess local pressure drop [ Δ P ( z ) P t ] in log scale, where we use the values of P t estimated from Fig. 11. The plots show a linear trend over one to three decades, depending on the capillary number, for the entire range of capillary numbers, with a linear range of the log –log plot that is all the larger as Ca is larger. The power-law exponents vary systematically with Ca within a range of 1–3 as shown by the blue and red straight lines in the plots. However, at very small values of the pressure drops, the data are noisy and deviate from the linearity, which corresponds to the gray symbols in the plots. This corresponds to the zone far behind the fingertip, where the Δ P ( z ) is very close to P t and the growth structure is almost frozen. As mentioned before, the statistical averaging becomes poorer as we move away from the fingertip, which therefore adds to the statistical noise in that region. We have therefore disregarded these points when fitting power laws to the data, and only rely on the points shown as colored symbols in the plots. The exponent β is thus measured from least-square fitting a power law to the colored data points. In Fig. 13, we plot the dependence of β on Ca for the two viscosity ratios, which exhibit a very similar trend, i.e., two plateaus for Ca 10 1.3 and Ca 10 0.2, and a crossover in between them. On the lower capillary number plateau, the values fluctuate within a range of 2.5–3.0. It then undergoes a crossover with the increase in Ca and then approaches another plateau toward 1. This is very similar to the rheological behavior of the widely studied steady-state flow, where the total two-phase flow rate varies non-linearly with the excess pressure drop at low capillary numbers, and then undergoes a crossover to a linear Darcy regime at high Ca.35–43 The origin of that non-linearity and the linear crossover has been well explained in terms of the pore-scale disorders in the capillary pressures,44,47 radii distribution,50 and wettabilities.51,52 Though the values of the steady-state non-linear exponent have been reported to differ for different systems, the range of them was similar to what we find here in Fig. 13.

FIG. 12.

Plot of the local excess pressure drop log 10 [ Δ ( P ( z ) P t ] between the invading and defending fluids as a function of the local growth rate of the finger, log 10 [ ν ( z ) ], for the viscosity ratios (a) M = 10 4 and (b) M = 10 5. Different sets correspond to different values of Ca as indicated by the symbols. The exponent β defined in Eq. (28) is calculated from the slopes. The slopes are calculated only for the range plotted with colored symbols, as the growth far behind the fingertip is almost frozen and leads to noisy data points as shown by the gray symbols. The blue and red straight lines are drawn with slopes 1 and 3.

FIG. 12.

Plot of the local excess pressure drop log 10 [ Δ ( P ( z ) P t ] between the invading and defending fluids as a function of the local growth rate of the finger, log 10 [ ν ( z ) ], for the viscosity ratios (a) M = 10 4 and (b) M = 10 5. Different sets correspond to different values of Ca as indicated by the symbols. The exponent β defined in Eq. (28) is calculated from the slopes. The slopes are calculated only for the range plotted with colored symbols, as the growth far behind the fingertip is almost frozen and leads to noisy data points as shown by the gray symbols. The blue and red straight lines are drawn with slopes 1 and 3.

Close modal
FIG. 13.

Dependence of the growth exponent β obtained from Fig. 12 on Ca. The error bars are of order 0.07 toward the lower side of Ca, whereas they are 0.01 on the higher side. They are too small compared to the scale of the plot, and therefore not shown. The two types of symbols correspond to the two values of M as indicated in the legend of the plot.

FIG. 13.

Dependence of the growth exponent β obtained from Fig. 12 on Ca. The error bars are of order 0.07 toward the lower side of Ca, whereas they are 0.01 on the higher side. They are too small compared to the scale of the plot, and therefore not shown. The two types of symbols correspond to the two values of M as indicated in the legend of the plot.

Close modal

We have addressed in this article how the growth rate of fingers during immiscible two-phase flow (drainage) in porous media depends on the local viscous pressure drops between the two fluids. In other words, we revisited the comparison of immiscible viscous fingers with Saffman–Taylor viscous fingering and DLA, for which a linear Laplacian growth is assumed. We considered pore networks with a regular geometry but with disorder in the pore/channel widths, and focused on the intermediate to higher capillary numbers, which include the intermediate regime between capillary and viscous fingers and the regime of pure viscous fingering. First, we looked into the geometrical properties of the drainage fingers and compared their statistical shape in the reference frame of the advancing finger to the smooth Saffman–Taylor fingers that develop in the continuum scale description of two-phase flow in porous media. Saffman and Taylor mentioned that viscous fingering in an empty Hele–Shaw cell without any porous structure is equivalent to what should be expected in a porous medium.12 Though this is true for single-phase flow, it is not necessarily true for two-phase flow, except perhaps if the porous medium is a regular lattice (no disorder), because the continuum description cannot account for the role of surface tension in a disordered porous medium. Indeed, in an empty Hele–Shaw cell, the fluid–fluid interface is continuous, and the long-range in-plane component of its curvature induces long-range forces, whereas in a porous medium, the interface is broken up into many small menisci, and surface tension acts at the scale of these menisci, which is the pore scale. The in-plane component of the curvature is then much larger, and, furthermore, stochasticity in the pore sizes plays an important role in shaping the displacement process by inducing a stochasticity in the capillary pressure thresholds.

We thus measured the statistical longitudinal volumetric density and growth rate profiles of the fractal drainage fingers in the dynamic reference frame attached to the most advanced fingertip, and compared them with the properties of smooth Saffman–Taylor fingers and of DLA growth structures. The width ratio with respect to the medium's width, λ, of the fingers' statistical density map, was found to approach the 0.5 Saffman–Taylor value at intermediate capillary numbers Ca, and then increase to around 0.8 when either increasing or decreasing the capillary numbers from that intermediate Ca range. The corresponding longitudinal volumetric density and growth rate profiles behave accordingly. For Saffman–Taylor fingers in continuum-medium, λ was observed to increase with the increase in surface tension,21,77,78 and it also increases with the angle of inlet wedge for systems with divergent wedges.29,55 In the case of our present system, it may be due to a combined effect of these two and needs further study.

The maximum growth happens at the most advanced fingertip and decreases exponentially behind it. Far behind the tip, the growth is almost frozen: the interfaces between the two fluids are held in place by the capillary barriers and thus do not move. We show that this growth behavior is directly correlated with the local pressure drop between the two phases, and controlled by the distribution of the capillary forces. By computing numerically an effective capillary threshold, we show that there exists a regime at intermediate capillary numbers where the linear Laplacian growth property does not hold as the local growth rate varies non-linearly with the excess local pressure drop across the interface (i.e., the excess between the local porous pressure drop and the capillary pressure threshold), Δ P. This non-linear regime of the growth of the invading finger is explained by accounting for the disorder distribution in the capillary thresholds in a theoretical assessment of the link between Δ P and the local growth rate. Indeed, when Δ P falls within the range of the capillary pressure thresholds, increasing the capillary number means that not only will the invasion velocity increase linearly with Δ P, but the number of pores along the front that are invaded at any time will also increase, thus rendering the growth rate non-linear. For a uniform capillary threshold distribution, the resulting growth law is expected to be quadratic. In our disorder pore networks, that distribution is not uniform, and consequently, the non-linear growth exponent β depends on the capillary number. The numerical simulations provide a range of β values similar to that of the rheology exponent previously reported for steady-state flow. In fact, it transitions from a plateau at values as large as 2.5–3 at the smallest investigated capillary numbers, to another plateau at 1 for the largest investigated capillary numbers. This large capillary number limit of β = 1 is to be expected, since it corresponds to configurations in which Δ P is larger than the upper boundary of the medium's capillary pressure thresholds; in such configurations, our theoretical assessment predicts a crossover of the flow regime to the linear Laplacian growth.

This study opens many prospects for future studies. One of them is the systematic investigation of the link between the capillary threshold distribution of the porous medium and the non-linear growth exponent β. Another one would be to study whether the values of β are comparable for steady-state flow and the present viscous fingering process in the same porous geometry. In order for the distribution of the capillary barriers to only depend on the pore-network geometry, the fluid–fluid interface must sample the entire barrier distribution at any time; otherwise, the capillary barrier/capillary threshold distribution will also depend on the flow patterns, in a manner that may differ between the steady-state flow and viscous fingering configurations.

This work is supported by the Research Council of Norway through its Centers of Excellence funding scheme, Project No. 262644.

The authors have no conflicts to disclose.

Santanu Sinha: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Investigation (equal); Methodology (equal); Project administration (equal); Software (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). Yves Méheust: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Project administration (equal); Writing – original draft (supporting); Writing – review & editing (equal). Hursanay Fyhn: Formal analysis (equal); Investigation (equal); Methodology (equal); Writing – review & editing (equal). Subhadeep Roy: Formal analysis (supporting); Investigation (equal); Methodology (equal); Writing – review & editing (equal). Alex Hansen: Conceptualization (equal); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Methodology (equal); Project administration (equal); Writing – original draft (supporting); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
R. L.
Chuoke
,
P.
van Meurs
, and
C.
van der Poel
, “
The instability of slow, immiscible, viscous liquid-liquid displacements in permeable media
,”
Trans. AIME
216
,
188
(
1959
).
2.
Y. C.
Yortsos
and
A. B.
Huang
, “
Linear-stability analysis of immiscible displacement. I. Simple basic flow profiles
,”
SPE Res. Eng.
1
,
378
(
1986
).
3.
D. W.
van Batenburg
,
S.
Berg
,
S.
Oedai
,
L. L.
David
,
A. O. N.
Siemens
, and
K.
Elewaut
, “
Visualisation of light oil mobilisation in ASP core floods using x-ray CT imaging
,”
paper presented at the SPE Asia Pacific Enhanced Oil Recovery Conference
,
2015
, Paper No. SPE-174660-MS.
4.
P.
Meakin
and
J. M.
Deutch
, “
The formation of surfaces by diffusion limited annihilation
,”
J. Chem. Phys.
85
,
2320
(
1986
).
5.
O. I.
Frette
,
K. J.
Måløy
,
J.
Schmittbuhl
, and
A.
Hansen
, “
Immiscible displacement of viscosity-matched fluids in two-dimensional porous media
,”
Phys. Rev. E
55
,
2969
(
1997
).
6.
T. K.
Perkins
,
O. C.
Johnston
, and
R. N.
Hoffman
, “
Mechanics of viscous fingering in miscible systems
,”
SPE J.
5
,
301
(
1965
).
7.
L.
Paterson
, “
Fingering with miscible fluids in a Hele Shaw cell
,”
Phys. Fluids
28
,
26
(
1985
).
8.
B.
Jha
,
L.
Cueto-Felgueroso
, and
R.
Juanes
, “
Fluid mixing from viscous fingering
,”
Phys. Rev. Lett.
106
,
194502
(
2011
).
9.
X.
Fu
,
L.
Cueto-Felgueroso
, and
R.
Juanes
, “
Viscous fingering with partially miscible fluids
,”
Phys. Rev. Fluids
2
,
104001
(
2017
).
10.
G.
Løvoll
,
Y.
Méheust
,
R.
Toussaint
,
J.
Schmittbuhl
, and
K. J.
Måløy
, “
Growth activity during fingering in a porous Hele-Shaw cell
,”
Phys. Rev. E
70
,
026301
(
2004
).
11.
I. C.
Salmo
,
K. S.
Sorbie
,
A.
Skauge
, and
M. A.
Alzaabi
, “
Immiscible viscous fingering: Modelling unstable water-oil displacement experiments in porous media
,”
Transp. Porous Media
145
,
291
(
2022
).
12.
P. G.
Saffman
and
G. I.
Taylor
, “
The penetration of a fluid into a porous medium or Hele-Shaw cell containing a more viscous liquid
,”
Proc. R. Soc. London
245
,
312
(
1958
).
13.
M. J.
Blunt
, “
Physically-based network modeling of multiphase flow in intermediate-wet porous media
,”
J. Pet. Sci. Eng.
20
,
117
(
1998
).
14.
J.
Feder
,
E. G.
Flekkøy
, and
A.
Hansen
,
Physics of Flow in Porous Media
(
Cambridge University Press
,
Cambridge
,
2022
).
15.
A.
Pinilla
,
M.
Asuaje
, and
N.
Ratkovich
, “
Experimental and computational advances on the study of viscous fingering: An umbrella review
,”
Heliyon
7
,
e07614
(
2021
).
16.
N.
Shahsavar
,
X.
Fu
, and
B.
Zhao
, “
Fingering instability during mixing-driven precipitation flow
,”
Transp. Porous Media
(published online
2023
).
17.
B.
Sandnes
,
E. G.
Flekkøy
,
H.
Knudsen
,
K. J.
Måløy
, and
H.
See
, “
Patterns and flow in frictional fluid dynamics
,”
Nat. Commun.
2
,
288
(
2011
).
18.
D.
Zhang
,
J. M.
Campbell
,
J. A.
Eriksen
,
E. G.
Flekkøy
,
K. J.
Måløy
,
C. W.
MacMinn
, and
B.
Sandnes
, “
Frictional fluid instabilities shaped by viscous forces
,”
Nat. Commun.
14
,
3044
(
2023
).
19.
F. K.
Eriksen
,
M.
Moura
,
M.
Jankov
,
A. L.
Turquet
, and
K. J.
Måløy
, “
Transition from viscous fingers to compact displacement during unstable drainage in porous media
,”
Phys. Rev. Fluids
7
,
013901
(
2022
).
20.
F.
Lanza
,
S.
Sinha
,
A.
Hansen
,
A.
Rosso
, and
L.
Talon
, “
Transition from viscous fingers to foam during drainage in heterogeneous porous media
,”
Phys. Fluids
35
,
103119
(
2023
).
21.
J. W.
McLean
and
P. G.
Saffman
, “
The effect of surface tension on the shape of fingers in a Hele Shaw cell
,”
J. Fluid Mech.
102
,
455
(
1981
).
22.
J.
Feder
,
Fractals
(
Springer
,
New York, NY
,
1988
).
23.
R.
Lenormand
,
E.
Touboul
, and
C.
Zarcone
, “
Numerical models and experiments on immiscible displacements in porous media
,”
J. Fluid Mech.
189
,
165
(
1988
).
24.
R.
Lenormand
and
C.
Zarcone
, “
Capillary fingering: Percolation and fractal dimension
,”
Transp. Porous Media
4
,
599
(
1989
).
25.
F.
Guo
and
S. A.
Aryana
, “
An experimental investigation of flow regimes in imbibition and drainage using a microfluidic platform
,”
Energies
12
,
1390
(
2019
).
26.
D.
Wilkinson
and
J. F.
Willemsen
, “
Invasion percolation: A new form of percolation theory
,”
J. Phys. A: Math. Gen.
16
,
3365
(
1983
).
27.
G. M.
Homsy
, “
Viscous fingering in porous media
,”
Annu. Rev. Fluid Mech.
19
,
271
(
1987
).
28.
K. J.
Måløy
,
J.
Feder
, and
T.
Jøssang
, “
Viscous fingering fractals in porous media
,”
Phys. Rev. Lett.
55
,
2688
(
1985
).
29.
A.
Arnéodo
,
Y.
Couder
,
G.
Grasseau
,
V.
Hakim
, and
M.
Rabaud
, “
Uncovering the analytical Saffman-Taylor finger in unstable viscous fingering and diffusion-limited aggregation
,”
Phys. Rev. Lett.
63
,
984
(
1989
).
30.
T. A. J.
Witten
and
L. M.
Sander
, “
Diffusion-limited aggregation, a kinetic critical phenomenon
,”
Phys. Rev. Lett.
47
,
1400
(
1981
).
31.
T. A. J.
Witten
and
L. M.
Sander
, “
Diffusion-limited aggregation
,”
Phys. Rev. B
27
,
5686
(
1983
).
32.
L.
Paterson
, “
Diffusion-limited aggregation and two-fluid displacements in porous media
,”
Phys. Rev. Lett.
52
,
1621
(
1984
).
33.
H.
Darcy
,
Les Fontaines publiques de la ville de Dijon
(
Victor Dalamont
,
Paris
,
1856
), p.
647
.
34.
S.
Whitaker
, “
Flow in porous media I: A theoretical derivation of Darcy's law
,”
Transp. Porous Media
1
,
3
(
1986
).
35.
K. T.
Tallakstad
,
H. A.
Knudsen
,
T.
Ramstad
,
G.
Løvoll
,
K. J.
Måløy
,
R.
Toussaint
, and
E. G.
Flekkøy
, “
Steady-state two-phase flow in porous media: Statistics and transport properties
,”
Phys. Rev. Lett.
102
,
074502
(
2009
).
36.
K. T.
Tallakstad
,
G.
Løvoll
,
H. A.
Knudsen
,
T.
Ramstad
,
E. G.
Flekkøy
, and
K. J.
Måløy
, “
Steady-state, simultaneous two-phase flow in porous media: An experimental study
,”
Phys. Rev. E
80
,
036308
(
2009
).
37.
O.
Aursjø
,
M.
Erpelding
,
K. T.
Tallakstad
,
E. G.
Flekkøy
,
A.
Hansen
, and
K. J.
Måløy
, “
Film flow dominated simultaneous flow of two viscous incompressible fluids through a porous medium
,”
Front. Phys.
2
,
63
(
2014
).
38.
E. M.
Rassi
,
S. L.
Codd
, and
J. D.
Seymour
, “
Nuclear magnetic resonance characterization of the stationary dynamics of partially saturated media during steady-state infiltration flow
,”
New J. Phys.
13
,
015007
(
2011
).
39.
E. M.
Rassi
,
S. L.
Codd
, and
J. D.
Seymour
, “
Corrigendum: Nuclear magnetic resonance characterization of the stationary dynamics of partially saturated media during steady-state infiltration flow
,”
New J. Phys.
16
,
039501
(
2014
).
40.
S.
Sinha
,
A. T.
Bender
,
M.
Danczyk
,
K.
Keepseagle
,
C. A.
Prather
,
J. M.
Bray
,
L. W.
Thrane
,
J. D.
Seymour
,
S. L.
Codd
, and
A.
Hansen
, “
Effective rheology of two-phase flow in three-dimensional porous media: Experiment and simulation
,”
Transp. Porous Media
119
,
77
(
2017
).
41.
Y.
Gao
,
Q.
Lin
,
B.
Bijeljic
, and
M. J.
Blunt
, “
Pore-scale dynamics and the multiphase Darcy law
,”
Phys. Rev. Fluids
5
,
013801
(
2020
).
42.
Y.
Zhang
,
B.
Bijeljic
,
Y.
Gao
,
Q.
Lin
, and
M. J.
Blunt
, “
Quantification of non-linear multiphase flow in porous media
,”
Geophys. Res. Lett.
48
,
e2020GL090477
, https://doi.org/10.1029/2020GL090477 (
2021
).
43.
Y.
Zhang
,
B.
Bijeljic
, and
M. J.
Blunt
, “
Nonlinear multiphase flow in hydrophobic porous media
,”
J. Fluid Mech.
934
,
R3
(
2022
).
44.
S.
Sinha
and
A.
Hansen
, “
Effective rheology of immiscible two-phase flow in porous media
,”
Europhys. Lett.
99
,
44004
(
2012
).
45.
S.
Sinha
,
A.
Hansen
,
D.
Bedeaux
, and
S.
Kjelstrup
, “
Effective rheology of bubbles moving in a capillary tube
,”
Phys. Rev. E
87
,
025001
(
2013
).
46.
H. L.
Cheon
,
H.
Fyhn
,
A.
Hansen
,
Ø.
Wilhelmsen
, and
S.
Sinha
, “
Steady-state two-phase flow of compressible and incompressible fluids in a capillary tube of varying radius
,”
Transp. Porous Media
147
,
15
(
2023
).
47.
S.
Roy
,
A.
Hansen
, and
S.
Sinha
, “
Effective rheology of two-phase flow in a capillary fiber bundle model
,”
Front. Phys.
7
,
92
(
2019
).
48.
A. G.
Yiotis
,
L.
Talon
, and
D.
Salin
, “
Blob population dynamics during immiscible two-phase flows in reconstructed porous media
,”
Phys. Rev. E
87
,
033001
(
2013
).
49.
S.
Roux
and
H. J.
Herrmann
, “
Disorder-induced nonlinear conductivity
,”
Europhys. Lett.
4
,
1227
(
1987
).
50.
S.
Roy
,
S.
Sinha
, and
A.
Hansen
, “
Role of pore-size distribution on effective rheology of two-phase flow in porous media
,”
Front. Water
3
,
709833
(
2021
).
51.
H.
Fyhn
,
S.
Sinha
,
S.
Roy
, and
A.
Hansen
, “
Rheology of immiscible two-phase flow in mixed wet porous media: Dynamic pore network model and capillary fiber bundle model results
,”
Transp. Porous Media
139
,
491
(
2021
).
52.
H.
Fyhn
,
S.
Sinha
, and
A.
Hansen
, “
Effective rheology of immiscible two-phase flow in porous media consisting of random mixtures of grains having two types of wetting properties
,”
Front. Phys.
11
,
1175426
(
2023
).
53.
R.
Toussaint
,
G.
Løvoll
,
Y.
Méheust
,
K. J.
Måløy
, and
J.
Schmittbuhl
, “
Influence of pore-scale disorder on viscous fingering during drainage
,”
Europhys. Lett.
71
,
583
(
2005
).
54.
A.
Arneodo
,
F.
Argoul
,
Y.
Couder
, and
M.
Rabaud
, “
Anisotropic Laplacian growths: From diffusion-limited aggregates to dendritic fractals
,”
Phys. Rev. Lett.
66
,
2332
(
1991
).
55.
A.
Arneodo
,
J.
Elezgaray
,
M.
Tabard
, and
F.
Tallet
, “
Statistical analysis of off-lattice diffusion-limited aggregates in channel and sector geometries
,”
Phys. Rev. E
53
,
6200
(
1996
).
56.
R.
Tao
,
M. A.
Novotny
, and
K.
Kaski
, “
Diffusion-limited aggregation with surface tension
,”
Phys. Rev. A
38
,
1019
(
1988
).
57.
M. A.
Novotny
,
R.
Tao
, and
D. P.
Landau
, “
Relaxation in DLA with surface tension
,”
J. Phys. A: Math. Gen.
23
,
3271
(
1990
).
58.
L.
Niemeyer
,
L.
Pietronero
, and
H. J.
Wiesmann
, “
Fractal dimension of dielectric breakdown
,”
Phys. Rev. Lett.
52
,
1033
(
1984
).
59.
C.
Amitrano
, “
Fractal dimensionality for the η model
,”
Phys. Rev. A
39
,
6618
(
1989
).
60.
J.
Mathiesen
and
M. H.
Jensen
, “
Tip splittings and phase transitions in the dielectric breakdown model: Mapping to the diffusion-limited aggregation model
,”
Phys. Rev. Lett.
88
,
235505
(
2002
).
61.
S.
Sinha
,
M. A.
Gjennestad
,
M.
Vassvik
, and
A.
Hansen
, “
Fluid meniscus algorithms for dynamic pore-network modeling of immiscible two-phase flow in porous media
,”
Front. Phys.
8
,
548497
(
2021
).
62.
E.
Aker
,
K. J.
Måløy
,
A.
Hansen
, and
G. G.
Batrouni
, “
A two-dimensional network simulator for two-phase flow in porous media
,”
Transp. Porous Media
32
,
163
(
1998
).
63.
M.
Erpelding
,
S.
Sinha
,
K. T.
Tallakstad
,
A.
Hansen
,
E. G.
Flekkøy
, and
K. J.
Måløy
, “
History independence of steady state in simultaneous two-phase flow through two-dimensional porous media
,”
Phys. Rev. E
88
,
053004
(
2013
).
64.
S.
Roy
,
S.
Sinha
, and
H.
Hansen
, “
Flow-area relations in immiscible two-phase flow in porous media
,”
Front. Phys.
8
,
4
(
2020
).
65.
S.
Roy
,
H.
Pedersen
,
S.
Sinha
, and
A.
Hansen
, “
The co-moving velocity in immiscible two-phase flow in porous media
,”
Transp. Porous Media
143
,
69
(
2022
).
66.
H.
Fyhn
,
S.
Sinha
, and
A.
Hansen
, “
Local statistics of immiscible and incompressible two-phase flow in porous media
,”
Physica A
616
,
128626
(
2023
).
67.
S.
Sinha
,
M. A.
Gjennestad
,
M.
Vassvik
,
M.
Winkler
,
A.
Hansen
, and
E. G.
Flekkøy
, “
Rheology of high-capillary number two-phase flow in porous media
,”
Front. Phys.
7
,
65
(
2019
).
68.
F. A. L.
Dullien
,
Porous Media: Fluid, Transport and Pore Structure
(
Academic Press
,
San Diego
,
1992
).
69.
E. W.
Washburn
, “
The dynamics of capillary flow
,”
Phys. Rev.
17
,
273
(
1921
).
70.
G. G.
Batrouni
and
A.
Hansen
, “
Fourier acceleration of iterative processes in disordered systems
,”
J. Stat. Phys.
52
,
747
(
1988
).
71.
D.
Wilkinson
, “
Percolation effects in immiscible displacement
,”
Phys. Rev. A
34
,
1380
(
1986
).
72.
V.
Joekar-Niasar
and
S. M.
Hassanizadeh
, “
Analysis of fundamentals of two-phase flow in porous media using dynamic pore-network models: A review
,”
Crit. Rev. Env. Sci. Technol.
42
,
1895
(
2012
).
73.
R.
Lenormand
and
C.
Zarcone
, “
Invasion percolation in an etched network: Measurement of a fractal dimension
,”
Phys. Rev. Lett.
54
,
2226
(
1985
).
74.
G.
Løvoll
,
Y.
Méheust
,
K. J.
Måløy
,
E.
Aker
, and
J.
Schmittbuhl
, “
Competition of gravity, capillary and viscous forces during drainage in a two-dimensional porous medium, a pore scale study
,”
Energy
30
,
861
(
2005
).
75.
E.
Somfai
and
R. C.
Ball
, “
Diffusion-limited aggregation in channel geometry
,”
Phys. Rev. E
68
,
020401(R)
(
2003
).
76.
Y.
Couder
,
F.
Argoul
,
A.
Arnéodo
,
J.
Maurer
, and
M.
Rabaud
, “
Statistical properties of fractal dendrites and anisotropic diffusion-limited aggregates
,”
Phys. Rev. A
42
,
3499
(
1990
).
77.
R.
Combescot
,
T.
Dombre
,
V.
Hakim
,
Y.
Pomeau
, and
A.
Pumir
, “
Shape selection of Saffman-Taylor fingers
,”
Phys. Rev. Lett.
56
,
2036
(
1986
).
78.
J.-M.
Vanden-Broeck
, “
Fingers in a Hele-Shaw cell with surface tension
,”
Phys. Fluids
26
,
2033
(
1983
).
79.
H.
Thomé
,
M.
Rabaud
,
V.
Hakim
, and
Y.
Couder
, “
The Saffman-Taylor instability: From the linear to the circular geometry
,”
Phys. Fluids
1
,
224
(
1989
).