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.

## I. INTRODUCTION

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 miscible^{6–9} and immiscible^{10,11} fluids in continuum^{1,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 compact^{19} or foam-like structures.^{20}

*zx*plane can be described by the following parametric equation given by Saffman and Taylor in their pioneering work in 1958:

^{12}

^{,}

*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:

*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 law

^{33,34}

*κ*and

*μ*are, respectively, the permeability of the porous medium and the fluid's viscosity. Taking into account the flow incompressibility, $ \u2207 \xb7 v = 0$ for incompressible fluids, the Laplace equation $ \u2207 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.

*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 $ \Delta P$ as a non-linear power law of exponent $ \u2248 1.85$.

^{35,36}Later, experiments with two incompressible fluids in the same porous medium found the exponent to be $ \u2248 1.35$ or 1.5 depending on the fractional flow.

^{37}Experiments with three-dimensional porous media made of glass bead packings

^{38–40}and real core samples,

^{41–43}performed by different groups, have further established this non-linearity with different exponent values. Various theoretical approaches

^{35,36,44,45}and numerical modelings with variable-radii tubes,

^{46}capillary bundles,

^{47}pore-networks models,

^{40,44}and lattice Boltzmann simulations

^{48}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*increase faster than the increase in $ \Delta P$, and hence, $ \beta > 1$.

^{49}When all the available pores along the interface, at any time, start flowing at a sufficiently high $ \Delta 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 sizes

^{50}and the wettabilities

^{51,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 media^{10,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 modeling^{61} 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.

## II. SYSTEM DESCRIPTION AND MODELING

^{61}which we tweaked to adapt to the present problem. The model has been developed for over a decade

^{62}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 \xd7 N L$ number of links embedded in a diamond lattice. Such a network with

*N*= 20 and

_{W}*N*= 80 is shown in Fig. 1. Here, the subscripts

_{L}*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 \u2009 mm$. The total network is therefore $ W = h N W \u2009 mm$ wide and $ L = h N L \u2009 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 \u2032$ as the interface (i.e., meniscus) between the two fluids moves along a link. We model this variation with a modified Young–Laplace equation

^{45,62,68}

*y*is the position of an interface inside a link

*i*. Such variation in the capillary pressure is shown in Fig. 2(b). Here, $ \gamma = \gamma \u2032 \u2009 cos \u2009 \theta $, where $ \gamma \u2032$ and

*θ*represent the surface tension between the two fluids and the contact angle, respectively. The

*r*is the average radius of the

_{i}*i*th 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.

^{23,69}

*i*th link. The network is assumed to be placed horizontally, and no gravity is considered. Here,

*μ*is the effective viscosity of the two fluids inside the link, which is given by $ \mu i = s i \mu n + ( 1 \u2212 s i ) \mu w$, where $ \mu n$ and $ \mu w$ are the viscosities of the non-wetting and wetting fluids, respectively, and

_{i}*s*is the non-wetting saturation in the

_{i}*i*th link, i.e., the proportion of the link volume that is occupied by the non-wetting fluid. The term

*k*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 = \pi r i 2$ is the cross-sectional area of the link. The summation of the capillary pressures $ p \u2032 i$ in Eq. (6) is over all the interfaces inside the link

_{i}*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 *p _{j}* at every node

*j*of the network, we use the Kirchhoff equation for incompressible flow $ \u2211 i \u2208 n j q i = 0$ for each node. Here the summation is over the links

*n*connected to a node

_{j}*j*. This provides a set of linear equations which we solve by conjugate-gradient method

^{70}with proper boundary conditions. This is illustrated in Fig. 1 where the two nodes (

*N*= 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

_{I}*Q*. We therefore set

*p*= 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.

_{j}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 \Delta t V j w / V j$ and $ V i n = q i \Delta 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 *q _{i}* in those links. Furthermore, technical details about this algorithm can be found in Ref. 61.

## III. RESULTS AND DISCUSSION

*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 as

^{12,71}

*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.

*Q*that must be imposed in the simulation for a given Ca, we first define it as

^{74}

*Q*

_{1}for the network saturated with wetting fluid, and then determine the flow rate

*Q*by using

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 *N _{W}* 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 $ \gamma = 0.03 \u2009 N / m$. Two viscosity ratios, $ M = \mu n / \mu w = 10 \u2212 4$ with $ \mu n = 10 \u2212 5 Pa s$ and $ \mu w = 10 \u2212 1 Pa s$, and $ M = 10 \u2212 5$ with $ \mu n = 10 \u2212 5 Pa s$ and $ \mu 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.

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 \u2212 4$ | 57 | 79 | 93 | 105 | 112 | 117 | 123 | 125 | 125 | 125 | 121 | 105 | 98 | 90 | 84 | 81 | 76 | 73 |

$ M = 10 \u2212 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 \u2212 4$ | 57 | 79 | 93 | 105 | 112 | 117 | 123 | 125 | 125 | 125 | 121 | 105 | 98 | 90 | 84 | 81 | 76 | 73 |

$ M = 10 \u2212 5$ | 30 | 49 | 56 | 66 | 71 | 74 | 81 | 81 | 82 | 81 | 80 | 65 | 62 | 57 | 54 | 51 | 46 | 45 |

### A. Longitudinal profiles of volume density and volumetric growth rate

*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*$ \sigma \u0303 ( x , z )$ as

*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

*h*

^{2}. 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 $ \sigma \u0303 ( x , z )$ map. In Fig. 3, we show the spatial distribution of $ \sigma \u0303 ( 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.

*growth rate*$ \nu \u0303 ( 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*),

*k*and

*k*– 1. Note that $ \nu \u0303 ( x , z )$ has the dimension of a velocity. By integrating $ \nu \u0303 ( x , z )$ in both

*x*and

*z*directions, we have

*Q*, we then define the

*growth density*as

*z*from the tip, by integrating $ \sigma \u0303 ( x , z ) , \u2009 \nu \u0303 ( x , z )$, and $ \varphi \u0303 ( x , z )$ in the

*x*direction,

*z*/

*h*, for different values of Ca and

*M*. Near the fingertips, the data show an exponential-type decay of $ \varphi ( z )$. This is indicated by the solid line where we plot a model function $ B \u2009 exp \u2009 ( \u2212 z h \xi )$ with $ \xi = 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 $ \u2212 1 / \xi $ and intercept $ ln \u2009 ( 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, $ \varphi ( 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.

*z*related to each other according to $ z = z \u0303 tip ( t ) \u2212 z \u0303 = v tip t \u2212 z \u0303$, where $ v tip$ is the velocity of the most advanced tip of the non-wetting fluid phase and $ z \u0303 tip ( t )$ is its position along $ z \u0303$. That velocity is constant in time, as shown in Fig. 5, where $ z \u0303 tip$ is plotted as a function of the time normalize by the time

*t*scaled by

*t*, where

_{f}*t*is the time at which the tip reaches the outlet boundary of the porous network. At any time

_{f}*t*, the volume of non-wetting fluid in a slice of the network normal to the mean flow direction and located between $ z \u0303$ and $ z \u0303 + \Delta z \u0303$ is $ \sigma \u0302 ( z \u0303 , t ) \Delta z \u0303$, with $ \Delta z \u0303 = \Delta z$. It results from the cumulative growth of the fingers between the time $ t 0 = z \u0303 / v tip$ at which the most advanced tip first arrived at $ z \u0303$ and

*t*, and can thus be expressed as

*x*,

*z*) reference frame, $ \sigma ( z )$, is related to $ \sigma \u0302 ( z \u0303 , t )$ through the following relation:

*z*direction as $ \Phi ( z ) = \u222b 0 z \varphi ( z \u2032 ) d z \u2032$ and $ \sigma c = Q / v tip$ leads to

To verify this functional form, the tip velocity is first measured from the slopes as the straight lines in Fig. 5 multiplied by *t _{f}*. We then plot $ \sigma ( z )$ scaled by

*σ*, for different values of Ca and

_{c}*M*(Fig. 6). Notice that, though $ \sigma ( z )$ is a function of the capillary numbers as shown in the insets, all $ \sigma ( z ) / \sigma 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 $ \Phi ( z )$, independent of Ca. We therefore calculate the average cumulative growth density $ \Phi ( 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 $ \sigma ( z ) / \sigma 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.

### B. The “continuum” shape

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 $ \sigma ( 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 $ \sigma \u0303 ( x , z ) \u2265 \sigma max / 2$ are shown. This region corresponds to the porous media with large invasion^{54,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 $ \lambda = 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.

*λ*, we measure a transverse volume density $ \rho ( 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

*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 DLA

^{55}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) $ \lambda = ( x + \u2212 x \u2212 ) / W$, where $ \rho ( x + ) = \rho ( x \u2212 ) = \lambda max / 2$ or alternatively from (B) $ \lambda = 1 / \rho 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 $ \u2248 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 $ \lambda \u2248 0.62$ for DLAs in linear channels

^{55,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 $ \lambda \u2248 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.

### C. Growth rate vs local pressure drop

*i*at (

*x*,

*z*) and

*q*(

*x*,

*z*) is the corresponding flow rate of the fluids in that link. Here, we consider $ \kappa = k \xaf / \mu w$, where $ k \xaf$ is the average mobility of the pores and

*μ*is the viscosity of the high-viscosity defending fluid. A pore will be invaded only when $ \Delta p ( x , z ) > P \u2032$, hence the Heaviside function $ \Theta ( \Delta p ( x , z ) \u2212 P \u2032 )$.

_{w}*z*only, we may average $ \Delta p ( x , z )$ over

*x*, which yields

*z*, $ [ \Delta P ( z ) \u2212 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 calculations

^{44}and observed by numerical simulations.

^{40,44,48}

^{50}and the pore wettabilities

^{51,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),

*β*is an exponent of non-linearity. As mentioned in the Introduction, such non-linear relationships with $ \beta > 1$ have been widely observed for steady-state flow in a certain regime of capillary numbers.

^{35–48}

*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

*z*, and

*P*

_{0}is the average wetting pressure at the most advanced fingertip,

*z*= 0. The pressure drop $ \Delta P ( z )$ between the wetting and non-wetting fluids at

*z*is then calculated as

In Fig. 10, $ P n ( z ) , \u2009 P w ( z )$, and $ \Delta 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 $ \Delta 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 $ \varphi ( z )$ in Fig. 4, which showed a rapid exponential fall within a similar distance, indicating that the growth $ \varphi ( z )$ is directly related to $ \Delta P ( z )$. Furthermore, we note that the plateau regimes of $ \Delta 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 $ \varphi ( 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 $ \Delta 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 $ \Delta P ( z )$ and then choose the one that provides the lowest least-square error for the straight-line fit for $ log \u2009 [ \nu ( z ) ]$ vs $ log \u2009 [ \Delta P ( z ) \u2212 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 $ \Delta 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.

Finally, we now set out to verify the agreement with Eq. (28) and determine the exponent *β*. In Fig. 12, we plot the growth rate $ \nu ( z )$ as a function of the excess local pressure drop $ [ \Delta P ( z ) \u2212 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 $ \Delta 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 \u2272 10 \u2212 1.3$ and $ Ca \u2273 10 \u2212 0.2$, and a crossover in between them. On the lower capillary number plateau, the values fluctuate within a range of $ \u2248 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.

## IV. DISCUSSION AND CONCLUSIONS

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), $ \Delta 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 $ \Delta P$ and the local growth rate. Indeed, when $ \Delta 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 $ \Delta 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 $ \Delta 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.

## ACKNOWLEDGMENTS

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

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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

## REFERENCES

*Physics of Flow in Porous Media*

*Les Fontaines publiques de la ville de Dijon*

*Porous Media: Fluid, Transport and Pore Structure*