Flow-Area Relations in Immiscible Two-Phase Flow in Porous Media

We present a theoretical framework for immiscible incompressible two-phase flow in homogeneous porous media that connects the distribution of local fluid velocities to the average seepage velocities. By dividing the pore area along a cut transversal to the average flow direction up into differential areas associated with the local flow velocities, we construct a distribution function that allows us to not only re-establish existing relationships of between the seepage velocities of the immiscible fluids, but also to find new relations between their higher moments. We support and demonstrate the formalism through numerical simulations using a dynamic pore-network model for immiscible two-phase flow with two- and three-dimensional pore networks. Our numerical results are in agreement with the theoretical considerations.

1. Introduction

When two immiscible fluids compete for the same pore space, we are dealing with immiscible two-phase flow in porous media [1]. A holy grail in porous media research is to find a proper description of immiscible two-phase flow at the continuum level, i.e., at scales where the porous medium may be treated as a continuum. Our understanding of immiscible two-phase flow at the pore level is increasing at a very high rate due to advances in experimental techniques combined with an explosive growth in computer power [2]. Still, the gap in scales between the physics at the pore level and a continuum description remains huge and the bridges that have been built so far across this gap are either complex to cross or rather rickety. To the latter class, we find the still dominating theory, first proposed by Wyckoff and Botset [3] and with an essential amendment by Leverett [4], namely relative permeability theory. The basic idea behind this theory is the following: Put yourself in the place of one of the two immiscible fluids. What does this fluid see? It sees a space in which it can flow limited by the solid matrix of the porous medium, but also by the other fluid. This reduces its mobility in the porous medium by a factor known as the relative permeability, which is a function of the how much space there is left for it. And here is the rickety part: this reduction of available space—expressed through the saturation—is the only parameter affecting the reduction factor or relative permeability. This is a very strong statement and clearly does not take into account that the distribution of immiscible fluid clusters will depend on how fast the fluids are flowing. Still, in the range of flow rates relevant for many industrial applications, this assumption works pretty well. It therefore, remains the essential work horse for practical applications.

Thermodynamically Constrained Averaging Theory (TCAT) [5–9] is built on the framework of relative permeability. However, it is based on a full analysis based on mechanical conservation laws, constitutive laws, e.g., for the motion of interfaces and contact lines, and on thermodynamics at the pore level. These are then scaled up using averaging theorems, which, loosely explained, consist of replacing derivatives of averages by averages of derivatives. In principle, this approach solves the up-scaling problem. However, as Gray and Miller point out in their book [9], each component of TCAT involves significant mathematical manipulations. The internal energy has contributions from the bulk liquids, the fluid-fluid and fluid-matrix interfaces at the pore level. The averaging process redefines the variables describing these contributions, but does not reduce their number. This accounts for a high level of complexity.

A further development somewhat along the same lines, based on non-equilibrium thermodynamics uses Euler homogeneity, more about this later, to define the up-scaled pressure. From this, Kjelstrup et al. derive constitutive equations for the flow [10, 11].

Another class of theories is based on detailed and specific assumptions concerning the physics involved. An example is Local Porosity Theory [12–17]. Another is DeProf theory which is a mechanical model combined with non-equilibrium statistical mechanics based on a classification scheme of fluid configurations at the pore level [18–20].

A recent work [21] has explored a new approach to immiscible two-phase flow in porous media based on elements borrowed from thermodynamics. That is, it is using the framework of thermodynamics, but without connecting it to processes involving heat. The spirit behind this approach is like that taken in Edwards and Oakshot's pseudo-thermodynamic theory of powders [22]. The approach consists in looking for general relations that transcends details of the physical processes involved. An example of such an approach in the field if immiscible two-phase flow in porous media is found in the Buckley-Leverett theory of invasion fronts [23]. The Buckley-Leverett equation is based solely on the principle of mass conservation and on the fractional flow rate being a function of the saturation. In the approach of Hansen et al.[21], equations are derived that originate from Euler homogeneity as in ordinary thermodynamics. These equations transcend the details of the physics involved in the same way that the equations of thermodynamics are universally applicable if a set of simple underlying conditions are met.

Thermodynamics is a theory that is valid on scales large enough so that the system it refers to may be regarded as a continuum. Statistical mechanics is then the theory that makes the connection between thermodynamics and the underlying atomistic picture.

It is the aim of this paper to formulate a description of immiscible two-phase flow in porous media that may form a link between the continuum-level approach of Hansen et al.[21] and the pore-level description of the problem—a sort of “statistical mechanics” from which the pseudo-thermodynamics may be derived, but which also describes the flow problem at the pore level.

After defining the system and the variables involved in section 2, we will in section 3 review the pseudo-thermodynamic approach [21]. The next section 4 we introduce the central object in the paper, the differential transversal area distribution which corresponds to the Boltzmann distribution in ordinary statistical mechanics, and relate it to the pseudo-thermodynamics relations. Then follows section 5 which then moves beyond the results of the pseudo-thermodynamics by focusing on fluctuations. In section 6 we use the dynamic network simulator [24] first introduced by Aker et al.[25] and then later refined [26–28] to verify the relations derived in the earlier sections. There is also a second goal behind this numerical work: the dynamic network model is a model at pore level and by its use, we show how the formalism developed here connect to the flow patterns at the pore level. Finally, we draw our conclusions in section 7.

2. System Definition

In two-phase flow, the steady state [29–31] is characterized by potentially strong fluctuations at the pore scale, but steady averages at the REV (Representative Elementary Volume) scale. As such they differ fundamentally from stationary states that are static at the pore scale as well. Steady states have much in common with ensembles in equilibrium statistical mechanics. They are also by implication assumed in the conventional descriptions of porous media flows that take the existence of an REV for granted.

Our REV is a block of homogeneous porous material of length L and area A. We prevent flow through the surfaces that are parallel to the L-direction which is the flow direction. The two remaining surfaces, each having an area A, act as inlet and outlet for the incompressible fluids that are injected and extracted from the REV. The porosity of the material is defined as

$\begin{array}{ll}\varphi =\frac{V_p}{V},& (1)\end{array}$

where V_p is the pore volume and V = AL is the volume of the REV. Due to the homogeneity of the porous medium, any cross section orthogonal to the axis along the L-direction will have a pore area that fluctuates around the value

$\begin{array}{ll}{A}_p=\frac{V_p}{L}=\varphi A.& (2)\end{array}$

There is also a solid matrix area fluctuating around

$\begin{array}{ll}{A}_s=A-A_p=(1-\varphi )A.& (3)\end{array}$

The homogeneity assumption consists in the fluctuations being so small that they can be ignored.

There is a time averaged volumetric flow rate Q through the REV. The volumetric flow rate consists of two components, Q_w and Q_n, which are the volumetric flow rates of the more wetting (w for “wetting”) and the less wetting (n for “non-wetting”) fluids with respect to the porous medium. They are related through

$\begin{array}{ll}Q=Q_w+Q_n.& (4)\end{array}$

In the porous medium, there is a volume V_w of the wetting fluid and a volume V_n of the non-wetting fluid so that V_p = V_w + V_n. We define the wetting and non-wetting saturations S_w = V_w/V_p and S_n = V_n/V_p, so that S_w + S_n = 1.

We define the wetting and non-wetting transversal pore areas A_w and A_n as the parts of the transversal pore area A_p which occupied by the wetting or the non-wetting fluids, respectively. We have that

$\begin{array}{ll}{A}_p=A_w+A_n.& (5)\end{array}$

As the porous medium is homogeneous, we will find the same averages A_w and A_n in any cross section through the porous medium orthogonal to the flow direction. We have therefore A_w/A_p = (A_wL)/(A_pL) = V_w/V_p = S_w, so that

$\begin{array}{ll}{A}_w=S_w{A}_p.& (6)\end{array}$

Likewise,

$\begin{array}{ll}{A}_n=S_n{A}_p=(1-S_w)A_p.& (7)\end{array}$

We define the seepage velocities, i.e., the average flow velocities in the pores, for the two immiscible fluids, v_w and v_n as

$\begin{array}{ll}{v}_w=\frac{Q_w}{A_w},& (8)\end{array}$

and

$\begin{array}{ll}{v}_n=\frac{Q_n}{A_n}.& (9)\end{array}$

The seepage velocity associated with the total flow rate Q is defined as

$\begin{array}{ll}{v}_p=\frac{Q}{A_p}.& (10)\end{array}$

We may express Equation (4) in terms of the seepage velocities,

$\begin{array}{ll}{v}_p=S_w{v}_w+S_n{v}_n.& (11)\end{array}$

3. Pseudo-Thermodynamic Relations

Hansen et al.[21] derived a number of relations between the seepage velocities defined in (8)–(10) based on the volumetric flow rate being an Euler homogeneous function of order one with respect to the wetting and non-wetting transversal pore areas A_w and A_n. We present here a short review of the main results in that paper for completeness. The meaning of the statement that the volumetric flow rate is an Euler homogeneous function of order one is that it obeys the scaling relation

$\begin{array}{ll}{Q}_p(\lambda A_w,\lambda A_n)=\lambda Q_p(A_w,A_n),& (12)\end{array}$

where λ is a scale factor. By taking the derivative of this equation with respect to λ and then setting λ = 1, we find

$\begin{array}{ll}{Q}_p(A_w,A_n)={\left(\frac{\partial Q_p}{\partial A_w}\right)}_{A_n}{A}_w+{\left(\frac{\partial Q_p}{\partial A_n}\right)}_{A_w}{A}_n.& (13)\end{array}$

By dividing this expression by the transversal pore area A_p and using Equations (5)–(7), we may write this equation as

$\begin{array}{ll}{v}_p=S_w{\left(\frac{\partial Q_p}{\partial A_w}\right)}_{A_n}+S_n{\left(\frac{\partial Q_p}{\partial A_n}\right)}_{A_w}.& (14)\end{array}$

The two partial derivatives have the units of velocity, and Hansen et al.[21] name these velocity functions the thermodynamic velocities,

$\begin{array}{ll}{\widehat{v}}_w={\left(\frac{\partial Q}{\partial A_w}\right)}_{A_n},& (15)\end{array}$

and

$\begin{array}{ll}{\widehat{v}}_n={\left(\frac{\partial Q}{\partial A_n}\right)}_{A_w}.& (16)\end{array}$

We use Equations (6) and (7) and the chain rule to derive

$\begin{array}{lll}\hfill {\left(\frac{\partial }{\partial A_w}\right)}_{A_n}& =& {\left(\frac{\partial S_w}{\partial A_w}\right)}_{A_n}{\left(\frac{\partial }{\partial S_w}\right)}_{A_p}\\ & +& {\left(\frac{\partial A_p}{\partial A_w}\right)}_{A_n}{\left(\frac{\partial }{\partial A_p}\right)}_{S_w}\\ & =& \frac{S_n}{A_p}{\left(\frac{\partial }{\partial S_w}\right)}_{A_p}+{\left(\frac{\partial }{\partial A_p}\right)}_{S_w}.& (17)\end{array}$

Likewise, we find

$\begin{array}{ll}{\left(\frac{\partial }{\partial A_n}\right)}_{A_w}=-\frac{S_w}{A_p}{\left(\frac{\partial }{\partial S_w}\right)}_{A_p}+{\left(\frac{\partial }{\partial A_p}\right)}_{S_w}.& (18)\end{array}$

We now combine these two equations with the definitions (15) and (16), and use that Q = A_pv_p, i.e., Equation (10), to find

$\begin{array}{ll}{\widehat{v}}_w=v_p+S_n\frac{d{v}_p}{d{S}_w},& (19)\end{array}$

and

$\begin{array}{ll}{\widehat{v}}_n=v_p-S_w\frac{d{v}_p}{d{S}_w}.& (20)\end{array}$

Combining the definitions (15) and (16) with Equation (14) gives

$\begin{array}{ll}{v}_p=S_w{\widehat{v}}_w+S_n{\widehat{v}}_n,& (21)\end{array}$

which should be compared to Equation (11). We see that

$\begin{array}{ll}{S}_w{v}_w+S_n{v}_n=S_w{\widehat{v}}_w+S_n{\widehat{v}}_n.& (22)\end{array}$

The seepage and thermodynamic velocities are related through a transformation $(v_w,v_n)\to ({\widehat{v}}_w,{\widehat{v}}_n)$ defining the co-moving velocity v_m,

$\begin{array}{ll}{\widehat{v}}_w=v_w+v_m{S}_n,& (23)\end{array}$

and

$\begin{array}{ll}{\widehat{v}}_n=v_n-v_m{S}_w.& (24)\end{array}$

We now calculate

$\begin{array}{lll}\hfill {\left(\frac{\partial Q}{\partial S_w}\right)}_{A_p}& =& {\left(\frac{\partial Q}{\partial A_w}\right)}_{A_n}{\left(\frac{\partial A_w}{\partial S_w}\right)}_{A_p}+\\ & & {\left(\frac{\partial Q}{\partial A_n}\right)}_{A_w}{\left(\frac{\partial A_n}{\partial S_w}\right)}_{A_p}.& (25)\end{array}$

Using Equations (6) and (7) together with Equations (19) and (20), we transform this equation into

$\begin{array}{ll}\frac{d{v}_p}{d{S}_w}={\widehat{v}}_w-{\widehat{v}}_n,& (26)\end{array}$

where we have used that v_p = Q/A_p, i.e., Equation (10). We now use Equation (21) to calculate

$\begin{array}{ll}\frac{d{v}_p}{d{S}_w}={\widehat{v}}_w-{\widehat{v}}_n+S_w\frac{d{\widehat{v}}_w}{d{S}_w}+S_n\frac{d{\widehat{v}}_n}{d{S}_w}.& (27)\end{array}$

Compare this equation to Equation (26) and we get an analog to the Gibbs-Duhem equation,

$\begin{array}{ll}{S}_w\frac{d{\widehat{v}}_w}{d{S}_w}+S_n\frac{d{\widehat{v}}_n}{d{S}_w}=0.& (28)\end{array}$

Using Equations (23) and (24), we find that the seepage velocities obey

$\begin{array}{ll}\frac{d{v}_p}{d{S}_w}=v_w-v_n+v_m,& (29)\end{array}$

and

$\begin{array}{ll}{S}_w\frac{d{v}_w}{d{S}_w}+S_n\frac{d{v}_n}{d{S}_w}=v_m,& (30)\end{array}$

where we have combined Equations (23) and (24) with Equation (28).

By combining Equations (15), (16), (23), and (24), one finds

$\begin{array}{ll}{v}_w=v_p+S_n\left(\frac{d{v}_p}{d{S}_w}-v_m\right),& (31)\end{array}$

and

$\begin{array}{ll}{v}_n=v_p-S_w\left(\frac{d{v}_p}{d{S}_w}-v_m\right).& (32)\end{array}$

These two equations, (31) and (32), may be seen as a transformation (v_p, v_m) → (v_w, v_n). The inverse of this transformation, i.e., (v_w, v_n) → (v_p, v_m) are given by Equations (11) and (29), i.e.,

$\begin{array}{lll}\hfill v_p& =& S_w{v}_w+S_n{v}_n,\\ \hfill v_m& =& S_w{v}_w^{\prime }+S_n{v}_n^{\prime },& (33)\end{array}$

where $v_w^{\prime }=d{v}_w/d{S}_w$ and $v_n^{\prime }=d{v}_n/d{S}_w$.

But, what is the co-moving velocity v_m physically? We first need to understand the thermodynamic velocities ${\widehat{v}}_w$ and ${\widehat{v}}_n$. These are the velocities the two fluids would have had if they were miscible. Equation (26) then tells us that a change in the saturation S_w leads to a change in the average seepage velocity v_p which is the difference in seepage velocities of the two fluids. However, the two fluids are not miscible and they do get in each other's way. How much is dictated by the co-moving velocity through Equation (29).

From Equation (26) onwards to the end of this sections, none of the equations contain the size of the REV. If we now imagine a REV associated with each point in the porous medium, we have a continuum description. We may then add equations that transport the fluids between these points. Assuming that the fluids are incompressible, these equations are [1]

$\begin{array}{ll}\varphi \frac{\partial S_w}{\partial t}=\frac{\partial \varphi S_w{v}_w}{\partial x},& (34)\end{array}$

where t is the time coordinate and x is the spatial coordinate, and

$\begin{array}{ll}\varphi \frac{\partial S_n}{\partial t}=\frac{\partial \varphi S_n{v}_n}{\partial x}.& (35)\end{array}$

We add the two equations and get

$\begin{array}{ll}\frac{\partial }{\partial x}\varphi v_p=0.& (36)\end{array}$

The generalization to three dimensions is straight forward.

In order to connect the equations that now have been derived to a given porous medium, constitutive equations for v_p and v_m need to be supplied, linking the flow to the driving forces. These may in the simplest case be pressure gradient and saturation gradient.

4. Differential Transversal Area Distributions

In this section, we connect the pseudo-thermodynamic results of section 3 to the properties of an underlying ensemble distribution. This concept in the context of immiscible two-phase flow was first considered by Savani et al.[32]. Here we generalize this concept. In some sense, we introduce here a statistical mechanics from which the pseudo-thermodynamics ensue.

We define a differential transversal pore area a_p = a_p(S_w, v) where v is a velocity such that a_pdv is the pore area covered by fluid, wetting or non-wetting, that has a velocity in the range [v, v + dv]. Hence, a_p—and the other differential transversal pore areas that we will proceed to construct—are statistical distributions of the pore level velocities. The new idea we are introducing is that the velocity distribution is measured in terms of transversal pore areas. This makes it possible to make the connection between the flow at the pore level and the pseudo-thermodynamic theory reviewed in the previous section.

We must have that

$\begin{array}{ll}{A}_p={\displaystyle {\int }_{-\infty }^{\infty }}dv{a}_p,& (37)\end{array}$

where the integral runs over the entire range of negative and positive velocities since there may be local areas where the flow direction is opposite to the global flow. The total flow rate Q is given by

$\begin{array}{ll}Q={\displaystyle {\int }_{-\infty }^{\infty }}dvv{a}_p,& (38)\end{array}$

and the see page velocity defined in Equation (10) is then given by

$\begin{array}{ll}{v}_p={\langle v\rangle }_p=\frac{1}{A_p}{\displaystyle {\int }_{-\infty }^{\infty }}dvv{a}_p.& (39)\end{array}$

Likewise, we define a wetting differential pore area a_w and a non-wetting differential pore area a_n. They have the same properties except that they are restricted to the wetting or the non-wetting fluids only. That is, we have

$\begin{array}{ll}{A}_w={\displaystyle {\int }_{-\infty }^{\infty }}dv{a}_w,& (40)\end{array}$

and

$\begin{array}{ll}{A}_n={\displaystyle {\int }_{-\infty }^{\infty }}dv{a}_n.& (41)\end{array}$

They relate to the wetting and non-wetting seepage velocities defined in Equations (8) and (9) as

$\begin{array}{ll}{v}_w={\langle v\rangle }_w=\frac{1}{A_w}{\displaystyle {\int }_{-\infty }^{\infty }}dvv{a}_w,& (42)\end{array}$

and

$\begin{array}{ll}{v}_n={\langle v\rangle }_n=\frac{1}{A_n}{\displaystyle {\int }_{-\infty }^{\infty }}dvv{a}_n.& (43)\end{array}$

We have that

$\begin{array}{ll}{a}_p=a_w+a_n.& (44)\end{array}$

We now combine this equation with Equation (39) to find

$\begin{array}{lll}\hfill v_p& =& \frac{1}{A_p}{\displaystyle {\int }_{-\infty }^{\infty }}dvv(a_w+a_n)\\ & =& \left(\frac{A_w}{A_p}\right)\frac{1}{A_w}{\displaystyle {\int }_{-\infty }^{\infty }}dvv{a}_w\\ & +& \left(\frac{A_n}{A_p}\right)\frac{1}{A_n}{\displaystyle {\int }_{-\infty }^{\infty }}dvv{a}_n\\ & =& S_w{v}_w+S_n{v}_n,& (45)\end{array}$

which is Equation (11). We have here used Equations (6) and (7).

We may associate a differential area a_m to the co-moving velocity v_m defined in Equation (29). By using Equations (39), (42), and (43) in combination with Equation (29), we find

$\begin{array}{lll}\hfill v_m& =& \frac{d{v}_p}{d{S}_w}-v_w+v_n\\ & =& \frac{1}{A_p}{\displaystyle {\int }_{-\infty }^{\infty }}dvv\left[\frac{\partial a_p}{\partial S_w}-\frac{a_w}{S_w}+\frac{a_n}{S_n}\right],& (46)\end{array}$

so that

$\begin{array}{lll}{a}_m& =& \frac{\partial a_p}{\partial S_w}-\frac{a_w}{S_w}+\frac{a_n}{S_n}\\ & =& \left(\frac{\partial a_w}{\partial S_w}-\frac{a_w}{S_w}\right)+\left(\frac{\partial a_n}{\partial S_w}+\frac{a_n}{S_n}\right),& (47)\end{array}$

where we have used Equation (44). Equation (47) may be rewritten as

$\begin{array}{ll}{a}_m=S_w\frac{\partial }{\partial S_w}\left(\frac{a_w}{S_w}\right)+S_n\frac{\partial }{\partial S_w}\left(\frac{a_n}{S_n}\right).& (48)\end{array}$

Averaging this equation over v and using Equations (42), (43), and (46) recovers Equation (30). Hence, we note that Equations (47) and (48) are the generalizations of Equations (29) and (30) to the differential transversal areas.

It follows that

$\begin{array}{ll}{A}_m={\displaystyle {\int }_{-\infty }^{\infty }}dv{a}_m=0,& (49)\end{array}$

where A_m is the pore area associated with co-moving velocity v_m. This is to expected as the areas A_w, A_n, A_p and A_m are ways to partition the transversal pore area A_p; and we have that A_w + A_n = A_p + 0. This implies that there is no volumetric flow rate associated with the co-moving velocity since

$\begin{array}{ll}{Q}_m=A_m{v}_m=0.& (50)\end{array}$

Lastly, we may associate differential transversal areas to the thermodynamic velocities defined in Equations (19) and (20). We use Equations (23) and (24) to find

$\begin{array}{ll}{\widehat{a}}_w=a_w+S_n{S}_w{a}_m,& (51)\end{array}$

and

$\begin{array}{ll}{\widehat{a}}_n=a_n-S_w{S}_n{a}_m,& (52)\end{array}$

where a_m is given in Equation (48). The thermodynamic velocities are then given by

$\begin{array}{ll}{\widehat{v}}_w=\frac{1}{A_w}{\displaystyle {\int }_{-\infty }^{\infty }}dvv{\widehat{a}}_w,& (53)\end{array}$

and

$\begin{array}{ll}{\widehat{v}}_n=\frac{1}{A_n}{\displaystyle {\int }_{-\infty }^{\infty }}dvv{\widehat{a}}_n.& (54)\end{array}$

We find as expected that

$\begin{array}{ll}{Â}_w={\displaystyle {\int }_{-\infty }^{\infty }}dv{\widehat{a}}_w=A_w,& (55)\end{array}$

and

$\begin{array}{ll}{Â}_n={\displaystyle {\int }_{-\infty }^{\infty }}dv{\widehat{a}}_n=A_n.& (56)\end{array}$

Summing the two differential transversal areas for the thermodynamic areas gives

$\begin{array}{ll}{\widehat{a}}_w+{\widehat{a}}_n=a_w+a_n=a_p.& (57)\end{array}$

This leads us to an important remark. The differential transversal areas are statistical velocity distributions at the pore level. We see that the differential transversal areas that are associated with the thermodynamic velocities are different from those associated with the seepage velocities. However, Equation (57) shows that the combined differential transversal area based upon the thermodynamic velocity distributions is the same as that based upon the distributions giving the seepage velocities. Hence, the two types of differential transversal areas represent a redistribution of the pore level velocities, but in such a way that A_w and A_n are preserved. We see the same from Equation (49) showing that A_m is zero and combining this Equations (51) and (52).

We see from Equation (47) that a_m is only zero if a_w and a_n are linear in S_w and S_n, respectively, i.e., a_w = S_wb_w where b_w is independent of S_w and a_n = S_nb_n where b_n is independent of S_n. Hence, this is the condition for the thermodynamic velocities to be equal to the seepage velocities.

5. Moments and Fluctuations

We define the qth moment of the seepage velocity distribution as

$\begin{array}{ll}{v}_p^q={\langle v^q\rangle }_p=\frac{1}{A_p}{\displaystyle {\int }_{-\infty }^{\infty }}dv{v}^q{a}_p.& (58)\end{array}$

By using Equation (44) we find immediately

$\begin{array}{ll}{v}_p^q=v_w^q{S}_w+v_n^q{S}_n,& (59)\end{array}$

where we have defined

$\begin{array}{ll}{v}_w^q={\langle v^q\rangle }_w=\frac{1}{A_w}{\displaystyle {\int }_{-\infty }^{\infty }}dv{v}^q{a}_w,& (60)\end{array}$

and

$\begin{array}{ll}{v}_n^q={\langle v^q\rangle }_n=\frac{1}{A_n}{\displaystyle {\int }_{-\infty }^{\infty }}dv{v}^q{a}_n.& (61)\end{array}$

We may work out the moments of the co-moving velocity are given by

$\begin{array}{ll}{v}_m^q=\frac{1}{A_p}{\displaystyle {\int }_{-\infty }^{\infty }}dv{v}^q{a}_m=\left[\frac{d{v}_p^q}{d{S}_w}-v_w^q+v_n^q\right],& (62)\end{array}$

where we have used (47).

The thermodynamic velocity moments may be defined as in a similar manner as the moments of the seepage velocities, (60) and (61),

$\begin{array}{ll}{\widehat{v}}_w^q={\langle {\widehat{v}}^q\rangle }_w=\frac{1}{A_w}{\displaystyle {\int }_{-\infty }^{\infty }}dv{v}^q{\widehat{a}}_w,& (63)\end{array}$

and

$\begin{array}{ll}{\widehat{v}}_n^q={\langle {\widehat{v}}^q\rangle }_n=\frac{1}{A_n}{\displaystyle {\int }_{-\infty }^{\infty }}dv{v}^q{\widehat{a}}_n.& (64)\end{array}$

and we find

$\begin{array}{ll}{\widehat{v}}_p^q={\widehat{v}}_w^q{S}_w+{\widehat{v}}_n^q{S}_n,& (65)\end{array}$

where we have used Equations (52) and (55).

We may Fourier transform a_p, a_w, and a_n,

$\begin{array}{ll}2\pi {\tilde{a}}_p(\omega )=A_p{\langle e^{iv\omega }\rangle }_p={\displaystyle {\int }_{-\infty }^{\infty }}dv{e}^{iv\omega }{a}_p,& (66)\end{array}$

$\begin{array}{ll}2\pi {\tilde{a}}_w(\omega )=A_w{\langle e^{iv\omega }\rangle }_w={\displaystyle {\int }_{-\infty }^{\infty }}dv{e}^{iv\omega }{a}_w,& (67)\end{array}$

and

$\begin{array}{ll}2\pi {\tilde{a}}_n(\omega )=A_n{\langle e^{iv\omega }\rangle }_n={\displaystyle {\int }_{-\infty }^{\infty }}dv{e}^{iv\omega }{a}_n.& (68)\end{array}$

From Equation (44) we find

$\begin{array}{ll}{\tilde{a}}_p(S_w,\omega )={\tilde{a}}_w(S_w,\omega )+{\tilde{a}}_n(S_w,\omega ),& (69)\end{array}$

and

$\begin{array}{ll}{\langle e^{iv\omega }\rangle }_p=S_w{\langle e^{iv\omega }\rangle }_w+S_n{\langle e^{iv\omega }\rangle }_n.& (70)\end{array}$

We write 〈exp(ivω)〉_p as a cumulant expansion,

$\begin{array}{ll}{\langle e^{iv\omega }\rangle }_p=\text{exp}\left({\displaystyle \sum _{k=1}^{\infty }}\frac{{(i\omega )}^k}{k!}{C}_p^k\right),& (71)\end{array}$

where $C_p^k$ is the kth cumulant. We define the wetting and non-wetting velocity cumulants $C_w^k$ and $C_n^k$ in the same way. We also write 〈exp(ivω)〉_p as a moment expansion

$\begin{array}{ll}{\langle e^{iv\omega }\rangle }_p={\displaystyle \sum _{m=0}^{\infty }}\frac{{(i\omega )}^m}{m!}{v}_p^m.& (72)\end{array}$

By expanding the cumulant expression in Equation (71) and equating each power in iω with the corresponding one in Equation (72), then repeating this for the wetting and non-wetting cumulants, and lastly combining them through Equation (70), we find for the term proportional to (iω)²,

$\begin{array}{ll}{C}_p^2+{(C_p^1)}^2=[C_w^2+{(C_w^1)}^2]S_w+[C_n^2+{(C_n^1)}^2]S_n.& (73)\end{array}$

Noting that $C_p^1=v_p$, $C_w^1=v_w$, and $C_n^1=v_n$ and using that $\Delta v_p^2=C_p^2$, $\Delta v_w^2=C_w^2$, and $\Delta v_n^2=C_n^2$, we find from this equation

$\begin{array}{ll}\Delta v_p^2=\Delta v_w^2{S}_w+\Delta v_n^2{S}_n+S_w{S}_n{(v_w-v_n)}^2.& (74)\end{array}$

We may follow this procedure for any of the cumulants.

The corresponding equation between the second cumulants of the thermodynamic velocities is

$\begin{array}{ll}\Delta {\widehat{v}}_p^2=\Delta {\widehat{v}}_w^2{S}_w+\Delta {\widehat{v}}_n^2{S}_n+S_w{S}_n{({\widehat{v}}_w-{\widehat{v}}_n)}^2.& (75)\end{array}$

6. Numerical Observations

The relations presented in sections 4, 5 provide the bridge between the velocity distributions at the pore level and the pseudo-thermodynamic theory outlined in section 3. In order to test these relations, and to show how they may be used, we use a dynamic pore network simulator [24].

In pore network modeling, the porous medium is represented by a network of pores which transport two immiscible fluids. The pore-network model we consider here can be applied to regular networks such as a regular lattice with an artificial disorder as well as to irregular networks such as a reconstructed network from real samples. The flow of the two immiscible fluids is described in this model by keeping the track of all interface positions with time. This approach of pore network modeling was first introduced by Aker et al.[25] for drainage displacements in a regular network. Over the last two decades, new mechanisms have been developed to extend the model for the steady-state flow as well as for irregular networks. A detailed description of this model in its most recent form can be found in Gjennestad et al. [26, 27] and Sinha et al. [28] and we therefore describe it here only briefly.

The porous medium is represented by a network of links that are connected at nodes. All the pore space in this model is assigned to the links and, hence, the nodes do not contain any volume, they only represent the positions where the links meet. The flow rate q_j inside any link j of the network at any instant of time for fully developed viscous flow is obtained by [33, 34],

$\begin{array}{ll}{q}_j=-\frac{g_j}{l_j{\mu }_j}{\displaystyle \left[\Delta p_j-\sum p_{c,j}\right]}& (76)\end{array}$

where Δp_j is the pressure drop across link, l_j is the link length and g_j is the link mobility which depends on the cross section of the link. The viscosity term μ_j is the saturation-weighted viscosity of the fluids inside the link given by μ_j = s_j,wμ_w + s_j,nμ_n where μ_w and μ_n are the wetting and non-wetting viscosities and s_j,w and s_j,n are the wetting and non-wetting fluid saturations inside the link, respectively. The term ∑p_c,j corresponds to the sum of all the interfacial pressures inside the jth link. A pore typically consists of two wider pore bodies connected by a narrow pore throat. We model this by using hour-glass shaped links. The variation of the interfacial pressure with the interface position for such a link is modeled by [34],

$\begin{array}{ll}|p_c(x)|=\frac{2\gamma \text{ cos }\theta }{r_j}\left[1-\text{cos}\left(\frac{2\pi x}{l_j}\right)\right]& (77)\end{array}$

where r_j is the average radius of the link and x ∈ [0, l_j] is the position of the interface inside the link. Here γ is the surface tension between the fluids and θ is the contact angle between the interface and the pore wall. These two Equations (77) and (76), together with the Kirchhoff relations, that is, the sum of the net volume flux at every node at each time step will be zero, provide a set of linear equations. We solve these equations with conjugate gradient solver [35] to calculate the local flow rates. All the interfaces are then advanced accordingly with small time steps. In order to achieve steady-state flow, we apply periodic boundary conditions in the direction of flow.

We construct a diamond lattice with 64 × 64 links in two dimensions (2D) with link lengths l_j = 1 mm for each link. Disorder is introduced by choosing the link radii r_j randomly from a uniform distribution in the range 0.1 mm and 0.4 mm. We use 10 different realizations of such network for our simulations in 2D. In three dimensions (3D), we use a network reconstructed from a 1.8 × 1.8 × 1.8 mm³ sample of Berea sandstone that contains 2, 274 links and 1, 163 nodes [36]. Simulations are performed under constant pressure drop ΔP across the network. For 2D, we have considered 3 different values for pressure drop such that, ΔP/L = 0.5, 1.0, and 1.5 MPa/m. For the 3D network, values of ΔP/L are chosen as, 10, 20, 40, and 80 MPa/m. The values for surface tension γ are chosen to be 0.02, 0.03, and 0.04N/m for both 2D and 3D. Three different values of viscosity ratios M(= μ_n/μ_w) = 0.5, 1.0, and 2.0 are considered. These values are chosen in such a way that the capillary number, defined as

$\begin{array}{ll}\text{Ca}=\frac{{\mu }_{e}Q}{\gamma A_p}& (78)\end{array}$

falls in a range of around 10⁻³ to 10⁻¹. Here μ_e is the saturation weighted effective viscosity of the system given by μ_e = S_wμ_w + S_nμ_n. Specifically, we find Ca in the range of 0.004–0.074 for 2D and 0.001–0.271 for 3D in the steady state. As the simulations are performed under constant pressure drop, the capillary number fluctuates. Ca is therefore calculated as functions of time by measuring the total flow rate Q along any cross section of the network perpendicular to the applied pressure drop. For any set of parameters, saturations are varied in the steps of 0.05 from 0 to 1 which correspond to 21 saturation values.

The simulations are continued to the steady state which is defined by the global measurable quantities, such as the fractional flow or the total flow rate Q fluctuate around a steady average. In the steady state, we calculate the seepage velocities averaged over time. First we use direct measurements, where we measure the global flow rates (Q, Q_w, and Q_n) and the pore areas (A_p, A_w, and A_n) through any cross section orthogonal to the applied pressure drop and then use Equations (8)–(10) to calculate the seepage velocities. Next, we perform the measurements using the differential pore areas (a_p, a_w, and a_n) and calculate the seepage velocities using description given in section 4. We then compare the results from the two measurements and calculate the co-moving velocities. We then verify the relation between the seepage velocities and their higher moments.

For the direct measurements, imagine a cross section at any place of the network orthogonal to the overall direction of flow. For the regular diamond lattice in 2D, all the links have the same length. Different moments of the seepage velocities can therefore be calculated by

$\begin{array}{ll}{v}_p^q=\frac{{\displaystyle \sum _j}{\left(\frac{q_j}{a_j}\right)}^q{a}_j}{{\displaystyle \sum _j}{a}_j},& (79)\end{array}$

$\begin{array}{ll}{v}_w^q=\frac{{\displaystyle \sum _j}{\left(\frac{q_j}{a_j}\right)}^q{a}_j{S}_{w,j}}{{\displaystyle \sum _j}{a}_j{S}_{w\text{,}j}},& (80)\end{array}$

and

$\begin{array}{ll}{v}_n^q=\frac{{\displaystyle \sum _j}{\left(\frac{q_j}{a_j}\right)}^q{a}_j{S}_{n\text{,}j}}{{\displaystyle \sum _j}{a}_j{S}_{n\text{,}j}},& (81)\end{array}$

where a_j is the projection of the pore area of the jth link on the cross sectional plane. Here, all links have the same angle α = 45° with the direction of the overall flow. However, in case of the irregular network in 3D, Equations (79)–(81) need to be modified as the links have different lengths and orientations. In such case, the one can calculate the seepage velocities by [28],

$\begin{array}{ll}{v}_p^q=\frac{{\displaystyle \sum _j}{\left(\frac{q_j}{a_j}\right)}^q{a}_j{l}_{x\text{,}j}}{{\displaystyle \sum _j}{a}_j{l}_{x\text{,}j}},& (82)\end{array}$

$\begin{array}{ll}{v}_w^q=\frac{{\displaystyle \sum _j}{\left(\frac{q_j}{a_j}\right)}^q{a}_j{S}_{w\text{,}j}{l}_{x\text{,}j}}{{\displaystyle \sum _j}{a}_j{S}_{w\text{,}j}{l}_{x\text{,}j}},& (83)\end{array}$

and

$\begin{array}{ll}{v}_n^q=\frac{{\displaystyle \sum _j}{\left(\frac{q_j}{a_j}\right)}^q{a}_j{S}_{n\text{,}j}{l}_{x\text{,}j}}{{\displaystyle \sum _j}{a}_j{S}_{n\text{,}j}{l}_{x\text{,}j}},& (84)\end{array}$

where l_x,j = l_j cos α_j is the projection of the link length (l_j) to the direction of the overall flow.

If we consider every link having the same length l_j = l and same orientations α_j = α in these equations, we retrieve the Equations (79)–(81). For the first moment (q = 1), the velocities $v_p^q$, $v_w^q$ and $v_n^q$ are equivalent to Q/A_p, Q/A_w, and Q/A_n, respectively, in both 2D and 3D.

For the second approach, we construct the distribution of differential transversal pore areas a_p, a_w, and a_n such that a_pdv, a_wdv, and a_ndv express the transversal pore areas for the total, wetting and non-wetting fluids within the velocity range from v to v + dv, so that they satisfy Equations (38), (40), and (41). We therefore have,

$\begin{array}{ll}\hfill a_p(v)dv& =\frac{1}{L}{\displaystyle \sum _j}{a}_j{l}_{x\text{,}j}\text{,}\\ \hfill a_w(v)dv& =\frac{1}{L}{\displaystyle \sum _j}{a}_j{l}_{x\text{,}j}{S}_{w\text{,}j}\text{,}\\ \hfill a_n(v)dv& =\frac{1}{L}{\displaystyle \sum _j}{a}_j{l}_{x\text{,}j}{S}_{n\text{,}j},& (85)\end{array}$

where j runs over all the sites satisfying the condition: v < v_j < v + dv, v_j being the local velocity of link j. In case of the 2D lattice, l_x,js are same for any j and given by $l_{x,j}=l/\sqrt{2}$. With these, different moments of the seepage velocities are then calculated using Equations (58), (60), and (61), respectively.

For any saturation, the seepage velocities and their higher moments should follow the relations (11), (45), and (59). We plot our numerical measurements in Figures 1, 2 for 2D and 3D, respectively. The upper row in each figure corresponds to the direct measurements and the lower row correspond to the measurements from the differential area distribution. A good agreement with the relations can be observed for first as well as for the higher moments for both the networks.

FIGURE 1

Figure 1. Verification of the relations (11), (45), and (59) between the steady-state seepage velocities v_p, v_w, and v_n, and their higher moments for the 2D regular network. The top row represents the direct approach of measurements using Equations (79), (80), and (81). The bottom row corresponds to the velocities measured from the differential area distributions defined in Equations (58), (60), and (61). v_p has a unit mm/s. Subsequently, for q = 2 and 3, the units for $v_p^q$ will be mm²/s and mm³/s, respectively.

FIGURE 2

Figure 2. Verification of the relations (11), (45), and (59) between the steady-state seepage velocities v_p, v_w, and v_n, and their higher moments for the 3D Berea network. The direct approach of measurements using Equations (82)–(84) are presented in the top row. The measurements using the differential area distributions defined in Equations (58), (60), and (61) are presented in the bottom row. v_p has a unit mm/s. Subsequently, for q = 2 and 3, the units for $v_p^q$ will be mm²/s and mm³/s, respectively.

Next we measure the fluctuations in the seepage velocities which obey Equation (74). Numerically, $\Delta v_p^2$, $\Delta v_w^2$, and $\Delta v_n^2$ are calculated from the knowledge of the 1st and 2nd moments by,

$\begin{array}{ll}\hfill \Delta v_p^2& ={\langle v^2\rangle }_p-{\langle v\rangle }_p^2,\\ \hfill \Delta v_w^2& ={\langle v^2\rangle }_w-{\langle v\rangle }_w^2,\\ \hfill \Delta v_n^2& ={\langle v^2\rangle }_n-{\langle v\rangle }_n^2.& (86)\end{array}$

In Figure 3, we plot these fluctuations for the two networks to compare with Equation (74) and good agreements is observed. There are some deviations in the results for the Berea network, since the results in 3D is based on only one network configuration whereas the results for 2D are averaged over 10 different configurations.

FIGURE 3

Figure 3. Numerical verification of Equation (74) between the fluctuations in the seepage velocities. $\Delta v_p^2$ has a unit mm²/s.

Finally, we verify the relations between seepage velocities and their higher moments while varying the fluid saturation as given by the Equations (29), (30), and (62). For this, we first calculated the co-moving velocity (v_m) and its higher moments from Equations (29) and (30) where we used the values of the seepage velocities measured with the direct approach. This is shown in the top rows of Figures 4, 5 for 2D and 3D, respectively, which show good agreements with Equations (29) and (30). We then compare these values of $v_m^q$ with the measurements from the differential transversal areas using Equation (46). For this, we first constructed the histogram for the differential pore area a_m corresponding to the co-moving velocity from Equation (47) where we have used the variations of a_p, a_w, and a_n with the saturation S_w. For this purpose, we have considered 21 different values of saturations within 0 and 1 with an interval of 0.05. We then integrate a_m from −∞ to ∞, weighted by the velocity and normalized by the total pore area to obtain the desired co-moving velocity with Equation (46). These results are plotted in the bottom row of Figures 4, 5 where they are compared with the results from direct measurements. The data points roughly follow the diagonal straight line showing satisfactory agreement with the theoretical formulations. However, we observe deviations in the results that is higher compared to the direct measurements. We believe this is due to the numerical errors that added up from several steps in the calculation such as the binning techniques while measuring the distributions, taking the derivatives and calculating the integrals. Moreover, the fluctuations for 3D are much higher compared to 2D, which is due to the lack of averaging over different samples as we have already mentioned earlier.

FIGURE 4

Figure 4. Measurement of the co-moving velocity (v_m) and its higher moments for the 2D network. The top row corresponds to the calculations using Equations (29) and (30) with the direct measurements. The bottom row shows the measurements of $v_m^q$ using the differential area distributions with Equation (46) and compared with the direct measurements where higher fluctuations are observed. v_m has a unit mm/s. Subsequently, for q = 2 and 3, the units for $v_m^q$ will be mm²/s and mm³/s, respectively.

FIGURE 5

Figure 5. Measurements of $v_m^q$ for the 3D Berea network where the top row corresponds to the direct measurements using Equations (29) and (30), and the bottom row corresponds to the measurement from the differential area distributions using Equation (46). Here, larger fluctuations in the results calculated with the differential pore area are observed compared to the 2D network. v_m has a unit mm/s. Subsequently, for q = 2 and 3, the units for $v_m^q$ will be mm²/s and mm³/s, respectively.

7. Summary

The aim of this paper is to provide the link between the pseudo-thermodynamic theory at the continuum level developed in Hansen et al. [21] (see section 3) and the velocities occurring at the pore level during immiscible two-phase flow in porous media. This link is provided by defining the differential transversal pore areas defined in section 4, which essentially correspond to the statistical distributions of velocities at the pore level. The central quantities are the velocity differential transversal pore area a_p, the wetting fluid differential velocity transversal pore area a_w, the non-wetting fluid velocity differential transversal pore area a_n, and the co-moving velocity differential transversal pore area a_m. We also consider the thermodynamic velocity differential transversal pore areas â_w and â_n. The relations found by Hansen et al.[21] for the average seepage velocities, the co-moving velocity and the thermodynamic velocities are generalized to the differential transversal areas here. In the following section 5, the relations are generalized to higher moments of the velocity distributions.

The theoretical derivations are then in section 6 validated by numerical simulations. We used dynamic pore-network modeling where an interface-tracking model is used to simulate steady-state two-phase flow. We used both regular pore networks and an irregular pore network reconstructed from a Berea sandstone for our simulations. By measuring the seepage velocities from the differential area distributions and comparing them with the direct measurements, we validated the essential predictions from the earlier theoretical sections.

Both Hansen et al.[21] and the present paper are to be seen as installments toward a theory for immiscible flow in porous media at the continuum scale. The structure of this theory reflects that found in thermodynamics: A set of general relations between the macroscopic variables based on energy conservation (i.e., the Gibbs relation) and Euler homogeneity. These general equations then have to be complemented by an equation of state which introduces the specifics of the system at hand. In the immiscible two-phase flow theory we are presenting here, Euler homogeneity and mass conservation provide the general equations that transcend the specifics of the porous medium. These general equations then have to be complemented by the constitutive equations for v_p and v_m, which provide the specifics of the porous medium.

The resulting set of equations may then be solved for structured porous media where the structure are associated with length scales larger than that set by the REV. This is e.g., seen in the explicit appearance of the porosity ϕ in Equations (34)–(36).

An open question, though, is what happens when there is non-trivial structure in the porous medium all the way from the pore scale to the continuum scale, see [37] and [38]—or when the saturation of the system is at a critical value, see [36]. The fundamental Euler scaling assumption (12) would then need to be modified, and with it, all the ensuing equations.

Data Availability Statement

The datasets generated for this study are available on request to the corresponding author.

Author Contributions

SR did the numerical simulations and analysis. SS developed the codes and performed the 3D simulations. AH developed the theory.

Funding

SS was supported by the National Natural Science Foundation of China under grant number 11750110430. This work was partly supported by the Research Council of Norway through its Center of Excellence funding scheme, project number 262644.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors thank Carl Fredrik Berg, M. Aa. Gjennestad, Knut Jørgen Måløy, Per Arne Slotte, Ole Torsæter, Morten Vassvik, and Mathias Winkler for interesting discussions. AH thanks the Beijing Computational Science Research Center CSRC for the hospitality and financial support.

References

1. Bear J. Dynamics of Fluids in Porous Media. Mineola, NY: Dover (1988).

Google Scholar

2. Blunt MJ. Multiphase Flow in Permeable Media. Cambridge: Cambridge University Press (2017).

3. Wyckoff RD, Botset HG. The flow of gas-liquid mixtures through unconsolidated sands. Physics. (1936) 7:325–45. doi: 10.1063/1.1745402

CrossRef Full Text | Google Scholar

4. Leverett MC. Capillary behavior in porous sands. Trans AIMME. (1940) 12:152.

Google Scholar

5. Hassanizadeh SM, Gray WG. Mechanics and thermodynamics of multiphase flow in porous media including interphase boundaries. Adv Water Res. (1990) 13:169–86.

Google Scholar

6. Hassanizadeh SM, Gray WG. Towards an improved description of the physics of two-phase flow. Adv Water Res. (1993) 16:53–67.

Google Scholar

7. Hassanizadeh SM, Gray WG. Thermodynamic basis of capillary pressure in porous media. Water Resour Res. (1993) 29:3389–405.

Google Scholar

8. Niessner J, Berg S, Hassanizadeh SM. Comparison of two-phase Darcy's law with a thermodynamically consistent approach. Transp Por Med. (2011) 88:133–48. doi: 10.1007/s11242-011-9730-0

CrossRef Full Text | Google Scholar

9. Gray WG, Miller CT. Introduction to the Thermodynamically Constrained Averaging Theory for Porous Medium Systems. Berlin: Springer Verlag (2014).

Google Scholar

10. Kjelstrup S, Bedeaux D, Hansen A, Hafskjold B, Galteland O. Non-isothermal transport of multi-phase fluids in porous media. The entropy production. Front Phys. (2018) 6:126. doi: 10.3389/fphy.2018.00126

CrossRef Full Text | Google Scholar

11. Kjelstrup S, Bedeaux D, Hansen A, Hafskjold B, Galteland O. Non-isothermal transport of multi-phase fluids in porous media. Constitutive equations. Front Phys. (2019) 6:150. doi: 10.3389/fphy.2018.00150

CrossRef Full Text | Google Scholar

12. Hilfer R, Besserer H. Macroscopic two-phase flow in porous media. Phys B. (2000) 279:125–9. doi: 10.1016/S0921-4526(99)00694-8

CrossRef Full Text | Google Scholar

13. Hilfer R. Capillary pressure, hysteresis and residual saturation in porous media. Phys A. (2006) 359:119–28. doi: 10.1016/j.physa.2005.05.086

CrossRef Full Text | Google Scholar

14. Hilfer R. Macroscopic capillarity and hysteresis for flow in porous media. Phys Rev E. (2006) 73:016307. doi: 10.1103/PhysRevE.73.016307

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Hilfer R. Macroscopic capillarity without a constitutive capillary pressure function. Phys A. (2006) 371:209–25. doi: 10.1016/j.physa.2006.04.051

CrossRef Full Text | Google Scholar

16. Hilfer R, Döster F. Percolation as a basic concept for capillarity. Transp Por Med. (2010) 82:507–19. doi: 10.1007/s11242-009-9395-0

CrossRef Full Text

17. Döster F, Hönig O, Hilfer R. Horizontal flow and capillarity-driven redistribution in porous media. Phys Rev E. (2012) 86:016317. doi: 10.1103/PhysRevE.86.016317

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Valavanides MS, Constantinides GN, Payatakes AC. Mechanistic model of steady-state two-phase flow in porous media based on Ganglion dynamics. Transp Porous Media. (1998) 30:267–99. doi: 10.1023/A:1006558121674

CrossRef Full Text | Google Scholar

19. Valavanides MS. Steady-state two-phase flow in porous media: review of progress in the development of the DeProF theory bridging pore-to statistical thermodynamics-scales. Oil Gas Sci Technol. (2012) 67:787–804. doi: 10.2516/ogst/2012056

CrossRef Full Text | Google Scholar

20. Valavanides MS. Review of steady-state two-phase flow in porous media: independent variables, universal energy efficiency map, critical flow conditions, effective characterization of flow and pore network. Transp Porous Media. (2018) 123:45–99. doi: 10.1007/s11242-018-1026-1

CrossRef Full Text | Google Scholar

21. Hansen A, Sinha S, Bedeaux D, Kjelstrup S, Gjennestad MA, Vassvik M. Relations between seepage velocities in immiscible, incompressible two-phase flow in porous media. Transp Porous Media. (2018) 125:565–87. doi: 10.1007/s11242-018-1139-6

CrossRef Full Text | Google Scholar

22. Edwards SF, Oakeshott RBS. Theory of powders. Phys A. (1989) 157, 1080–90.

Google Scholar

23. Buckley SE, Leverett MC. Mechanism of fluid displacements in sands. Trans AIME. (1942) 146:107–17.

Google Scholar

24. Joekar-Niasar V, Hassanizadeh SM. Analysis of fundamentals of two-phase flow in porous media using dynamic pore-network models: a review. Crit Rev Environ Sci Technol. (2012) 42: 1895–76. doi: 10.1080/10643389.2011.574101

CrossRef Full Text | Google Scholar

25. Aker E, Måløy KJ, Hansen A, Batrouni GG. A two-dimensional network simulator for two-phase flow in porous media. Transp Porous Media. (1998) 32:163–86. doi: 10.1023/A:1006510106194

CrossRef Full Text | Google Scholar

26. Gjennestad MA, Vassvik M, Kjelstrup S, Hansen A. Stable and efficient time integration of a dynamic pore network model for two-phase flow in porous media. Front Phys. (2018) 6:56. doi: 10.3389/fphy.2018.00056

CrossRef Full Text | Google Scholar

27. Gjennestad MA, Winkler M, Hansen A. Pore network modeling of the effects of viscosity ratio and pressure gradient on steady-state incompressible two-phase flow in porous media. arXiv:1911.07490 (2019).

Google Scholar

28. Sinha S, Gjennestad MA, Vassvik M, Hansen A. A dynamic network simulator for immiscible two-phase flow in porous media. arXiv:1907.12842 (2019).

Google Scholar

29. Tallakstad KT, Knudsen HA, Ramstad T, Løvoll G, Måløy KJ, Toussaint R, et al. Steady-state two-phase flow in porous media: statistics and transport properties. Phys Rev Lett. (2009) 102:074502. doi: 10.1103/PhysRevLett.102.074502

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Tallakstad KT, Løvoll G, Knudsen HA, Ramstad T, Flekkøy EG, Måløy KJ. Steady-state simultaneous two-phase flow in porous media: an experimental study. Phys Rev E. (2009) 80:036308. doi: 10.1103/PhysRevE.80.036308

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Aursjø O, Erpelding M, Tallakstad KT, Flekkøy EG, Hansen A, Måløy KJ. Film flow dominated simultaneous flow of two viscous incompressible fluids through a porous medium. Front Phys. (2014) 2:63. doi: 10.3389/fphy.2014.00063

CrossRef Full Text | Google Scholar

32. Savani I, Bedeaux D, Kjelstrup S, Sinha S, Vassvik M, Hansen A. Ensemble distribution for immiscible two-phase flow in porous media. Phys Rev E. (2017) 95:023116. doi: 10.1103/PhysRevE.95.023116

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Washburn EW. The dynamics of capillary flow. Phys Rev. (1921) 17:273. doi: 10.1103/PhysRev.17.273

CrossRef Full Text | Google Scholar

34. Sinha S, Hansen A, Bedeaux D, Kjelstrup S. Effective rheology of bubbles moving in a capillary tube. Phys Rev E.(2013) 87:025001. doi: 10.1103/PhysRevE.87.025001

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Batrouni GG, Hansen A. Fourier acceleration of iterative processes in disordered systems. J Stat Phys. (1988) 52:747–73. doi: 10.1007/BF01019728

CrossRef Full Text | Google Scholar

36. Ramstad T, Hansen A, Øren PE. Flux-dependent percolation transition in immiscible two-phase flow in porous media. Phys Rev E. (2009) 79:036310. doi: 10.1103/PhysRevE.79.036310

CrossRef Full Text | Google Scholar

37. Parteli EJR, da Silva LR, Andrade JS Jr. Self-organized percolation in multi-layered structures. J Stat Mech. (2010) 2010:P03026. doi: 10.1088/1742-5468/2010/03/P03026

CrossRef Full Text | Google Scholar

38. Hansen A, da Silva LR, Lucena L. Spatial correlations in permeability distributions due to extreme dynamics restructuring of unconsolidated sandstone. Phys A. (2011) 390:553. doi: 10.1016/j.physa.2010.10.011

CrossRef Full Text | Google Scholar