Relations between Seepage Velocities in Immiscible, Incompressible Two-Phase Flow in Porous Media

Based on thermodynamic considerations we derive a set of equations relating the seepage velocities of the fluid components in immiscible and incompressible two-phase flow in porous media. They necessitate the introduction of a new velocity function, the co-moving velocity. This velocity function is a characteristic of the porous medium. Together with a constitutive relation between the velocities and the driving forces, such as the pressure gradient, these equations form a closed set. We solve four versions of the capillary tube model analytically using this theory. We test the theory numerically on a network model.


Introduction
The simultaneous flow of immiscible fluids through porous media has been studied for a long time [1]. It is a problem that lies at the heart of many important geophysical and industrial processes. Often, the length scales in the problem span numerous decades; from the pores measured in micrometers to reservoir scales measured in kilometers. At the largest scales, the porous medium is treated as a continuum governed by effective equations that encode the physics at the pore scale.
The problem of tying the pore scale physics together with the effective description at large scale is the upscaling problem. In 1936, Wycoff and Botset proposed a generalization of the Darcy equation to immiscible two-phase flow [2]. It is instructive to reread Wycoff and Botset's article. This is where the concept of relative permeability is introduced. The paper is eighty years old and yet it is still remarkably modern. Capillary pressure was first considered by Richards as early as 1931 [3]. In 1940, Leverett combined capillary pressure with the concept of relative permeability, and the framework which dominates all later practical analyses of immiscible multiphase flow in porous media was in place [4].
It is the aim of this paper to present a new theory for flow in porous media that is based on thermodynamic considerations. In the same way as Buckley and Leverett's analysis based on the conservation of the mass of the fluids in the porous medium led to their Buckley-Leverett equation [23], the thermodynamical considerations presented here leads to a set of equations that are general in that they transcend detailed assumptions about the flow. This is in contrast to the relative permeability equations that rely on a number of specific physical assumptions.
The theory we present focuses on one aspect of thermodynamic theory. We see it as a starting point for a more general analysis based on non-equilibrium thermodynamics [24,25,26] which combines conservation laws with the laws of thermodynamics. The structure of the theory, already as it is presented here, is reminiscent of the structure of thermodynamics itself: we have a number of variables that are related through general thermodynamic principles leaving an equation of state to account for the detailed physics of the problem. To our knowledge, the analysis we present here has no predecessor. However, the framework in which it resides is that originally laid out by Gray and Hassanizadeh [6], Hassanizadeh and Gray [7,8,9], and Gray and Miller [12].
The theory we present concerns relations between the flow rates, or equivalently, the seepage velocities of the immiscible fluids. We do not discuss relations between the seepage velocities and the driving forces that create them, such as the pressure gradient. In this sense, we are presenting a "kinetic" theory of immiscible two-phase flow in porous media, leaving out the "dynamic" aspects which will be addressed elsewhere. Our main result is a set of equations between the seepage velocities that together with a constitutive relation between the average seepage velocity and the driving forces lead to a closed set of equations.
We consider only one-dimensional flow in this paper, deferring to later the generalization to three dimensions. The fluids are assumed to be incompressible.
In section 2 we describe the porous medium system we consider. We review the key concepts that will be used in the subsequent discussion. In particular, we discuss the relation between average seepage velocity and the seepage velocities of each of the two fluids.
In section 3, we introduce our thermodynamic considerations. We focus on the observation that the total volumetric flow rate is an Euler homogeneous function of the first order. This allows us to define two thermodynamic velocities as derivatives of the total volumetric flow. We then go on in section 4 to deriving several equations between the thermodynamic velocities, one of which is closely related to the Gibbs-Duhem equation in thermodynamics. In section 4.1, we point out that the seepage velocity of each fluid is generally not equal to the corresponding thermodynamic velocity. This is due to the constraints that the geometry of the porous medium puts on how the immiscible fluids arrange themselves. In ordinary thermodynamics, such constraints are not present and questions of this type do not arise. We relate the thermodynamic and seepage velocities through the introduction of a co-moving velocity function. This co-moving velocity is a characteristic of the porous medium and depends of the driving forces only through the velocities and the saturation.
We discuss in section 5 what happens when the thermodynamic and seepage velocities coincide, when the wetting and non-wetting seepage velocities coincide and when the thermodynamic wetting and non-wetting velocities coincide. Normally, the coincidences appear at different saturations. However, under certain conditions, the three coincide. When this happens, the fluids behave as if they were miscible.
In section 6 we write down the full set of equations to describe immiscible two-phase flow in porous media. These equations are valid on a large scale where the porous medium may be seen as a continuum.
In section 7, we analyze four versions of the capillary tube model [27,28] within the concept of the representative elementary volume (REV) developed in the previous sections. This allows us to demonstrate these concepts in detail and to demonstrate the internal consistency of the theory.
In section 8 we analyze a network model [29] for immiscible two-phase flow within the framework of our theory. We calculate the co-moving velocity for the model using a square and a hexagonal lattice. The co-moving velocity is to within the level of the statistical fluctuations equal for the two lattice topologies. We compare successfully the measured seepage velocities for each fluid component with those calculated from our theory. Lastly, we use our theory to predict the coincidence of the three saturations defined in section 5, and verifying this prediction numerically. Section 9 summarizes the results from the previous sections together with some further remarks on the difference between our approach and that of relative permeability. We may nevertheless anticipate here what will be our conclusions. Through our theory we have accomplished two things. The first one is to construct a closed set of equations for the seepage velocities as a function of saturation based solely on thermodynamic principles. We do not propose relations between seepage velocities and driving forces such as pressure gradient here. Our discussion is solely based on relations between the seepage velocities. The second accomplishment is to pave the way for further thermodynamic analysis by identifying the proper thermodynamic variables that relate to the flow rates.

Defining the system
The aim of this paper is to derive a set of equations on the continuum level where differentiation make sense. We define a representative elementary volume -REV -as a block of porous material with no internal structure filled with two immiscible and incompressible fluids: it is described by a small set of parameters which we will now proceed to define.
We illustrate the REV in Fig. 1. It is a block of homogeneous porous material of length L and area A. We seal off the surfaces that are parallel to the L-direction (the flow direction). The two remaining surfaces, each with area A, are kept open and act as inlet and outlet for the fluids that are injected and extracted from the REV. The porosity is where V p is the pore volume and V = AL the total volume of the REV. Due to the homogeneity of the porous medium, any cross section orthogonal to the axis along the L-direction (named the x axis for later) will reveal a pore area that fluctuates around the value while the solid matrix area fluctuates around 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 (see figure 1). 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. We have All flows have as frame of reference the non-deformable solid matrix. In the porous medium, there is a volume V w of incompressible wetting fluid and a volume V n of incompressible 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 pore areas A w and A n as the parts of the pore area A p which are filled with the wetting or the non-wetting liquids respectively. 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  Fig. 1 In the upper part of the figure, we see the REV from the side. There is a flow Q = Qw + Qn across it. An imaginary cut is made through the REV in the direction orthogonal to the flow. In the lower left corner, the surface of the imaginary cut is illustrated. A magnification of the surface of the cut is shown in the lower right corner. The pore structure is illustrated as brown and black circles. The pores that are brown, are filled with wetting fluid and the pores that are black, are filled with non-wetting fluid. The wetting fluid-filled pores form in total an area Aw and the non-wetting fluid-filled pores form in total an area An. The total pore area of the imaginary cut in the lower left corner is Ap = Aw + An = Aφ.
orthogonal to the flow direction. This is illustrated in Fig. 1. We have that Likewise, where We define the seepage velocities for the two immiscible fluids, v w and v n as and v n = Q n A n .
Hence, equation (4) may be written We finally define a seepage velocity associated with the total flow rate Q as By using equations (6) to (8) and (9), (10) and (12) we transform (11) into All the variables in this equation can be measured experimentally. The extensive variables describing the REV, among them the wetting and non-wetting pore areas A w and A n , and the wetting and non-wetting volumetric flow rates Q w and Q n , are averages over the REV and their corresponding intensive variables S w , S n , v w and v n are therefore by definition a property of the REV.
This definition of a REV is the same as that used by Gray and Miller [12] in their thermodynamically constrained averaging theory and in the earlier literature on thermodynamics of multi-phase flow, see e.g. [7]. The REV can be regarded as homogeneous only on the REV scale. With one driving force the variables on this scale can be obtained, knowing the micro-scale ensemble distribution [30,31]. This averaging procedure must keep invariant the entropy production, a necessary condition listed already in the 1990ies [7,12].

The volumetric flow rate Q is an Euler homogeneous function of order one
The volumetric flow rate Q is a homogeneous function of order one of the areas A w and A n defined in Section 2. We now proceed to demonstrate this. Suppose we scale the three areas A w → λA w , A n → λA n and A s → λA s , keeping in mind that which we get from equations (2) and (3). This corresponds to enlarging the area A of the REV to λA as illustrated in figure 2. This transformation does not change the porosity φ. Due to equation (14) with φ kept fixed, the only two independent variables in the problem are then A w and A n , and as a consequence the volumetric flow rate depends on these two variables and not on A s , Q = Q(A w , A n ). The scaling A w → λA w and A n → λA n affects the volumetric flow rate as follows: As long as the porous material does not have a fractal structure, this scaling property is essentially self evident as explained in Fig. 2.
We take the derivative with respect to λ on both sides of (15) and set λ = 1. This gives 1

Equation (16) is essentially the Euler theorem for homogeneous functions of order one. By dividing this equation by
where we have used equations (6) and (7). The two partial derivatives in equation (17) have the units of velocity. Hence, we define two thermodynamic velocitiesv w andv n aŝ Equation (17) The fact that the volumetric flow rate Q is a homogeneous function of order one in the variables A w and A n implies that the thermodynamic velocities  (18) and (19) must be homogeneous of the zeroth order so that havev This implies that they must depend on the areas A w and A n through their ratio A w /A n = S w /S n = S w /(1 − S w ), where we have used equations (6) -(8) and we may writev

Consequences of the Euler theorem: new equations
We now change the variables from We divide by the area A p and find replacing the partial derivative by a total derivative in accordance with equation (22). We then proceed by taking the derivative of equation (20) with respect to S w , Combining this equation with equation (24), we find This equation is a Gibbs-Duhem-like equation, here relating the thermodynamic velocitiesv w andv n . We note that equations (20), (24) and (26) are interrelated in that any pair selected from the three equations will contain the third.
The thermodynamic velocities (18) and (19) may be expressed in terms of the variables A p and S w , where we have used that Likewise, we findv We introduce a function g = g(S w ) defined as Combining this equation with equation (26) gives We now take the derivative of equation (24) with respect to S w and combine with the two previous equations (30) and (31) to find where we have used equation (5). Hence, we have finally and

New equations in terms of the seepage velocities
Equation (20) has the same form as equation (13), This does not imply thatv w = v w andv n = v n . The most general relations we can write down betweenv w and v w , andv n and v n and still fulfill (35) arê andv where v m is a function of S w with the units of velocity. We name this function the co-moving velocity. We may interpret the physical contents of the second equality in (35) in the following way: The thermodynamic velocitiesv w and v n differ from the v w and v n as each fluid carries along some of the other fluid -hence "co-moving." To preserve the volume, these co-moving contributions are related in such a way that A n S w v m is the co-moving volume flow of the non-wetting fluid with the wetting fluid and A w S n v m is the counter-moving volume flow of the wetting fluid with the non-wetting fluid. By using equations (36) and (37), we find that the seepage velocities (9) and (10) also are homogeneous of the zeroth order so that Equations (20), (24) and (26) are in terms of the thermodynamic velocities defined by the area derivatives of the volumetric flow, (18) and (19). We now express them in terms of the seepage velocities v w and v n defined in equations (9) and (10). We do this by invoking the transformations (36) and (37). Equation (24) then becomes where we have used equations (36) and (37) and v m is the co-moving velocity. Expressing (26) in terms of the seepage velocities (9) and (10), we find The three equations (13), (39) and (40) are dependent: combine any two of them and the third will follow.
By a derivation similar to that which led to equations (33) and (34), we find and

Cross points
We may define three wetting saturations of special interest, S A w , S B w and S C w , that have the following properties: (36) and (37), we have that v m = 0 for this wetting saturation. We furthermore have that v m ≤ 0 for S w ≤ S A w and v m ≥ 0 for S w ≥ S A w .
-B When S w = S B w , we have thatv w =v n . From equation (35) we find v w =v n = v. Equations (27) and (29) make dv/dS w = 0 equivalent to the equality between the three velocities at this wetting saturation. Furthermore, at dv/dS w = 0, v has its minimum value. We have that dv/dS w ≤ 0 for S w ≤ S B w and dv/dS w ≥ 0 for S w ≥ S B w .
-C When S w = S C w , we have that v w = v n . This implies that v w = v n = v from equation (13) for this wetting saturation. By equations (45) and (46) we have that this is equivalent to v m = dv/dS w .
By considering the relations between the velocity variables derived so far, we may show that we either have or It is straight forward to demonstrate that if either two of the three wetting saturations S A w , S B w or S C w , coincide, then the third must also have the same value. If this is the case, we have v w = v n and dv/dS w = 0, so that the two immiscible fluids behave as if they were miscible [38]. We will use our theory to identify the wetting saturation at which this coincidence occurs for the network model we study in section 8 and then verify numerically that this is indeed so.

A closed set of equations
We now consider the porous medium on scales larger than the REV scale. The properties of the porous medium may at these scales vary in space. We consider here only flow in the x direction, deferring the generalization to higher dimensions to later. Let x be a point somewhere in this porous medium. Hence, all the variables in the following will be functions of x.
Equations (13) and (40), where v ′ w = dv w /dS w and v ′ n = dv n /dS w , can be seen as a transformation from the velocity pair (v w , v n ) to the pair (v, v m ). We invert the velocity transformation, and where v ′ = dv/dS w . We now turn to the equations controlling the time evolution of the flow. The transport equation for S w is Likewise, we have for S n Summing these two equations and using equations (5) and (13), we find the incompressibility condition ∂v ∂x = 0 .
Equations (13) In order to close the equation set, we need an equation for the co-moving velocity v m . This is a constitutive equation that characterizes how the immiscible fluids affect each other due to the constraints imposed by the porous Specifying v m will specify the form of the differential equation (40). We will work out v m analytically for four systems in section 7 and map it for a range of parameters in a network model in section 8. Note that it is not enough to specify S w and v to determine v m . It is also a function of dv/dS w , since dv/dS w depends on how the external parameters are controlled as S w is changed. For example, if the total volumetric flow rate Q = A p v is held constant when S w is changed, then dv/dS w = 0 for all values of S w , and the system traces the curve (S w , v, v ′ = 0) in (S w , v, v ′ ) space. If, on the other hand, the driving forces are held constant when S w is changed, then v and dv/dS w will follow some non-trivial curve in (S w , v, v ′ ) space, see section 8. A third example is to control the wetting and the non-wetting volumetric flow rates Q w and Q n , making S w and v and dv/dS w dependent variables, following curves in (S w , v, v ′ ) space depending on how Q w and Q n are changed.
The driving forces enter the discussion through the constitutive equation where f represents the driving forces, e.g., f = −∂P/∂x, where dP/dx is the pressure gradient. The equation set (13), (40), (47) and (49) together with the constitutive equation for v m is independent of the driving forces. They are valid for any constitutive equations (51) that may be proposed. In particular, non-linear constitutive equations such as those recently proposed which suggest that two immiscible fluids behave as if they were a single Bingham plastic [32,33,34,35,36] fits well into this framework. This is in stark contrast to e.g. the relative permeability framework that assumes a particular constitutive equation for v w and v n containing the two relative permeabilities k r,w (S w ) and k r,n (S w ) in addition to the capillary pressure function P c (S w ). The equation set we present here is thus much more general than that of the relative permeability framework.

Analytically tractable models
We will in this section analyze four variants of the capillary tube model [27,28]. In each case, we calculate the co-moving velocity v m and then use it to demonstrate the consistency of our theory. We also use this section to clarify the physical meaning of the area derivatives introduced in section 3.

Parallel capillaries filled with either fluid
In this simplest case, we envision a bundle of N parallel capillaries, all equal. Each tube has a cross-sectional inner (i.e. pore) area a p . N w of these capillaries are filled with the wetting fluid and N n are filled with the non-wetting fluid, so that N w + N n = N . Let us assume that the seepage velocity in the capillaries filled with wetting fluid is v w and the seepage velocity in the capillaries filled with non-wetting fluid is v n . The total volumetric flow rate is then The wetting saturation is given by We calculate the derivative of Q with respect to A w to determine the thermodynamic velocityv w defined in equation (18) together with equation (52), The derivative is thus performed changing the number of capillaries from N to N + δN and letting all the added capillaries contain the wetting fluid so that N w → N w + δN w = N w + δN . Likewise, we find using equation (19) v n = ∂Q ∂A n Aw = 1 a p δN n [Q(N w , N n ) − Q(N w , N n + δN n )] = v n . (55) From equations (36) and (37) we find that the co-moving velocity is zero, v m = 0 (56) for this system. Hence, we have that the seepage velocities v w and v n are equal to the thermodynamic velocitiesv w andv n .

Parallel capillaries with bubbles
Suppose that we have N parallel capillaries as shown in Fig. 3. Each tube has a length L and an average inner area a p . The diameter of each capillary varies along the long axis. Each capillary is filled with a bubble train of wetting or non-wetting fluid. The capillary forces due to the interfaces vary as the bubble train moves due to the varying diameter. We furthermore assume that the wetting fluid, or more precisely, the more wetting fluid (in contrast to the non-wetting fluid which is the less wetting fluid ) does not form films along the pore walls so that the fluids do not pass each other, see [37] for details. The volume of the wetting fluid in each tube is L w a p and the volume of the nonwetting fluid is L n a p . Hence, the saturations are S w = L w /L and S n = L n /L for each tube. Suppose now that the seepage velocity in each tube is v when averaged over time. Both the wetting and non-wetting seepage velocities must be equal to the average seepage velocity since the bubbles do not pass each other: We now make an imaginary cut through the capillaries orthogonal to the flow direction as shown in Fig. 3. There will be a number N w capillaries where the cut passes through the wetting fluid and a number N n capillaries where the cut passes through the non-wetting fluid. Averaging over time, we must have and Hence, the wetting and non-wetting areas defined in section 2 are and A n = N n a p ; .
We also have that We will now calculate the derivative (18) defining the thermodynamic velocityv w . We do this by changing the pore area A p → A p +δA p = N a p +a p δN . We wish to change A w while keeping A n fixed. This can only be done by adjusting S w while changing δN . This leads to the two equations and We solve either (63) or (64) for δS w (they contain the same information), finding We recognize that this equation is equation (27). Hence, we could have taken this equation as the starting point of discussing this model. However, doing the derivative explicitly demonstrates its operational meaning.
Rather than performing a similar calculation forv n as in (66) forv w , we use equation (29) and havê We now combine equations (57), (66) and (67) with (36) and (37) and read We see that the expression for v m is independent of the constitutive equation that relates flow rate to the driving forces. Rather, it expresses the relation between the thermodynamic velocitiesv w andv n , and the seepage velocities, reflecting the co-moving -or rather the lack of it -in the capillaries.
This result gives us an opportunity to clarify the physical meaning of the co-moving velocity v m . We have assumed that the bubbles cannot pass each other as shown in Fig. 3. Equation (68) gives upon substitution in equations (36) and (37), and v n =v n + S w dv dS w .
We see that these two equations compensate exactly for the two equations (27) and (29),v that give the thermodynamic velocities resulting in v w = v n = v.

Parallel capillaries with a subset of smaller ones
We now turn to the third example. There are N parallel capillaries of length L. A fraction of these capillaries has an inner area a s , whereas the rest has an inner area a l . The total pore area carried by the small capillaries is A s and the total pore area carried by the larger capillaries is A l so that We define the fraction of small capillaries to be S w,i , so that We now assume that the small capillaries are so narrow that only the wetting fluid enters them. These pores constitute the irreducible wetting-fluid contents of the model in that we cannot go below this saturation. However, they still contribute to the flow. Hence, at a saturation S w which we assume to be larger than or equal to S w,i which is then the irreducible wetting fluid saturation, we have that the wetting pore area is where we have defined We assume there is a seepage velocity v sw in the small capillaries, a seepage velocity v lw in the larger capillaries filled with wetting fluid and a seepage velocity v n in the larger filled with non-wetting fluid. The velocities v sw , v lw and v n are independent of each other. The total volumetric flow rate is then given by or in term of average seepage velocity We will now calculate the derivative (18) defining the thermodynamic velocityv w . We could have done this using equation (27). However, it is instructive to perform the derivative yet again through area differentials so that their meaning become clear.
In order to calculate the derivative (19) we change the pore area A p → A p + δA p so that δA w = δA p and δA n = 0. In addition, we have that δA s = S w,i δA p , so that δA w = δA s + δA lw = S w,i δA p + δA lw = δA p .
Hence, we have We now combine this expression with equation (74) for the total volumetric flow Q to find We could have foundv n in the same way. However, using (29) combined with (75) givesv We now use equation (37) to find As a check, we now use the co-moving velocity found in (81) together with equation (36) which is the expected result.

Large capillaries with bubbles and small capillaries with wetting fluid only
The fourth example that we consider is a combination of the two previous models discussed in section 7.2 and 7.3, namely a set of N capillary tubes. A fraction S w,i of these capillaries has an inner area a s , whereas the rest has an inner area a l . The capillaries with the larger area contain bubbles as in section 7.2. The smaller capillaries contain only wetting liquid as in section 7.3. The fluid contents in the small capillaries is irreducible. The seepage velocity in the smaller capillaries is v sw . The wetting and non-wetting seepage velocities in the larger capillaries are the same, Hence, the seepage velocity is then (Ca) Fig. 4 The co-moving velocity vm may be expressed as in equation (39): vm = v ′ − vw + vn, or as in equation (40): vm = Swv ′ w + Snv ′ n . We plot here the two different expressions for vm against each other. The data points are shaded according to the capillary number Ca of the flow. The data are based both on the hexagonal and the square lattice.
We have that dv since v sl is independent of S w . Equation (29) giveŝ where we have used (84). We now use equations (37) and (83) Combining equations (85), (86) and (87) gives Hence, the co-moving velocity is a sum of the co-moving velocities of the two previous sections, see (68) and (81).
We check the consistency of this calculation by using the co-moving velocity found in (88) together with equation (45) to calculate v w . We find as expected These four analytically tractable examples have allowed us to demonstrate in detail how the thermodynamic formalism that we have introduced in this paper works. We have in particular calculated the co-moving velocity v m in all four cases. As is clear from these examples, it does not depend on the constitutive equation or any other equation except through v and dv/dS w .
Hence, the discussion of these four models has not included the constitutive equation (51). In order to have the seepage velocities v w and v n as functions of the driving forces, we would need to include the constitutive equation and the incompressibility condition (49) in the analysis.

Network model studies
Our aim in this section is to map the co-moving velocity v m for a network model which we describe in the following. We compare the measured seepage velocities with those calculated using the formulas derived earlier. We use the theory to predict where the three cross points defined and discussed in section 5 coincide and verify this prediction through direct numerical calculation.
The network model first proposed by Aker et al. [29] has been refined over the years and is today a versatile model for immiscible two-phase flow under steady-state conditions due to the implementation of bi-periodic boundary conditions. The model tracks the interfaces between the immiscible fluids by solving the Kirchhoff equations with a capillary pressure in links created by the interfaces they contain due to a surface tension σ. The links are hour glass shaped with average radii r drawn from a probability distribution. The minimum size of the bubbles of wetting or non-wetting fluid in a given link has a minimum length of r, the average radius of that link.
Our network and flow parameters were chosen as follows: we used hexagonal (honeycomb) lattices of size 60 × 40 and square lattices of size 80 × 30, both with all links having the same length l = 1 mm. The radii r of the links were assigned from a flat distribution in the interval 0.1l < r < 0.4l. The surface tension σ was either 0.02, 0.03 or 0.04 N/m. The viscosities of the fluids η w and η n were set to either either 10 −3 , 2 × 10 −3 , 5 × 10 −3 and 10 −2 Pa s. The pressure drop over the width of the network, |∆P/L| was held constant during each run and set to either 0.1, 0.15, 0.2 0.3, 0.5 or 1 M Pa/m. The capillary number Ca ranged from around 10 −3 to 10 −1 , see figure 4.
We measure the seepage velocities v w as follows. In both networks, we define layers normal to the flow direction. Each layer has an index j and is the set of links intersected by a plane normal to the flow direction. In layer j, each link has an index i. For every time step, we calculate where S w,ij is the saturation of link i in layer j and a w,ij is the cross-sectional area of the link, projected into the plane normal to the flow direction. The flow rate q w,ij is measured as the volume of wetting fluid that has flowed past We show in figure 4 (40)). This verifies the main results of section 4.1, besides being a measure of the quality of the numerical derivatives based on central difference for the internal points. Forward/backward difference was used at the endpoints of the curves where S w = 0 or S w = 1. Figure 5 shows the co-moving velocity v m calculated from equation (40) as a function of S w and v ′ for the flow parameters given above. The data are roughly consistent with a linear form where a ≈ −0.15, b ≈ 0.79 and c ≈ −0.095. We note that we show both the data from the hexagonal and the square lattices in figure 5. The co-moving velocity v m appears not to be sensitive to the lattice topology. We compare in figure 6 the measured seepage velocities v w and v n to the seepage velocities calculated from equations (45) and (46) where we have used the co-moving velocity shown in figure 5. The spikes at low v w and v n are due to the transition from both fluids moving to only one fluid moving creating a jump in the derivative.
We note that when the wetting saturation S w is equal to S w,0 = c/a in equation (91), the co-moving velocity v m takes the form v m = bv ′ . For the flows we study here, we find S w,0 ≈ 0.6. With reference to the discussion of the cross points in section 5, we then have v m = 0 and v ′ = 0 occurring for the same wetting saturation. Hence, we have S A w = S B w = S C w at this point. Figure 5 shows two closely spaced curves, one dotted and one dashed. The dashed curve passes through the point (v m , v ′ ) = (0, 0) (marked by a non-filled circle). We show in figure 7 v, v w , v n ,v w andv n corresponding to these two curves. The thermodynamic velocities have been calculated using equations (27) and (29). In 7a, the flow parameters have been chosen so that v ′ = 0 and v m = 0 do not coincide as we see in figure 5. The ordering of the cross points then follow the inequality (44): S C w ≤ S B w ≤ S A w . In figure 7b, we have chosen the parameters so that v ′ = v m = 0 for a given saturation, S w = S w,0 ≈ 0.6, see figure 5. The theory then predicts S C w = S B w = S A w , which is exactly what we find numerically. This constitutes a non-trivial test of the theory: from an experimental point of view, this is the point where v w = v n = v and v has its minimum. These are all measurable quantities.

Discussion and conclusion
We have here introduced a new thermodynamic formalism to describe immiscible two-phase flow in porous media. This work is to be seen as a first step towards a complete theory based on equilibrium and non-equilibrium thermodynamics. In fact, the only aspect of thermodynamic analysis that we have utilized here is the recognition that the central variables are Euler homogeneous functions. This has allowed us to derive a closed set of equations in section 6. This was achieved by introducing the co-moving velocity (50). This is a velocity function that relates the seepage velocities v w and v n to the thermodynamic velocitiesv w andv n defined through the Euler theorem in section 3.
The theory we have developed here is "kinetic" in the sense that it only concerns relations between the velocities of the fluids in the porous medium and not relations between the driving forces that generate these velocities and the velocities. We may use the term "dynamic" about these latter relations.
The theory rests on a description on two levels, or scales. The first one is the pore scale. Then there is REV scale. This scale is much larger than the pore scale so that the porous medium is differentiable. At the same time, we  figure  5 passing through (vm, v ′ ) = (0, 0). In this case we find that the cross points all coincide. Hence, we have S C w = S B w = S A w , as predicted.
assume the REV to be so small that we may assume that the flow through it is homogeneous and at steady state. At scales much larger than the REV scale, we do not assume steady-state flow. Rather, saturations and velocity fields may vary both in space and time. It is at this scale that the equations in section 6 apply. The necessity to introduce the co-moving velocity v m is surprising at first. However, it is a result of the homogenization procedure when going from the pore scale to the REV scale: The pore geometry imposes restrictions on the flow which must reflect itself at the REV scale.
One may at this point ask what has been gained in comparison to the relative permeability formalism? The relative permeability approach consists in making explicit assumptions about the functional form of v w and v n through the generalized Darcy equations [1]. The relative permeability formalism thus reduces the immiscible two-phase problem to the knowledge of three functions k r,w = k r,w (S w ), k r,n = k r,w (S n ) and P c = P c (S w ). However, strong assumptions about the flow have been made in order to reduce the problem to these three functions; assumptions that are known to be at best approximative.
Our approach leading to the central equations in section 6, does not make any assumptions about the flow field beyond the scaling assumption in equation (15). Hence, the constitutive equation (51) can have any form, for example the non-linear one presented in [32,33,34,35,36]. Hence, our approach is in this sense much more general than the relative permeability one.