Flow of immiscible fluids in porous media [1–4] is a problem that has been in the hands of engineers for a long time. This has resulted in a schism between the physics at the pore scale and the description of flow at scales where the porous medium may be seen as a continuum, also known as the Darcy scale. This fact has led to the phenomenological relative permeability equations proposed by Wyckoff and Botset in 1936 [5] with the inclusion by Leverett of the concept of capillary pressure in 1940 [6], still being the unchallenged approach to calculate flow in porous media at large scales. The basic ideas of the theory are easy to grasp. Seen from the viewpoint of one of the immiscible fluids, the solid matrix and the other fluid together reduce the pore space in which that fluid can move. Hence, the effective permeability as seen by the fluid is reduced, and the permeability reduction factor of each fluid is their relative permeability. The capillary pressure models the interfacial tension at the interfaces between the immiscible fluids by assuming that there is a pressure difference between the pressure fields in each fluid. The central variables of the theory, relative permeability, and capillary pressure curves are determined routinely in the laboratory under the name special core analysis and then used as input in reservoir simulators [3, 7]. A key assumption in the theory is that the relative permeability and the capillary pressure are functions of the saturation alone. This simplifies the theory tremendously, but it also distances the theory from realism.

Some progress has been made in order to improve on relative permeability theory. Barenblatt et al. [8] recognized that a key assumption in relative permeability theory is that locally, the fluids are in phase equilibrium, even if the flow as a whole is developing. This has as an implication that the central variables of that theory, relative permeability, and the capillary pressure are functions of the saturation alone. They then go on to generalize the theory to flow which is locally out of equilibrium, exploring how the central variables change. Wang et al. [9] went further by introducing dynamic length scales due to the mixing zone variations, over which the spatial averaging is performed.

The relative permeability approach is phenomenological and going beyond relative permeability theory means making a connection between the physics at the pore scale with the Darcy scale description of the flow. The attempts at constructing a connection between these two scales—which may stretch from micrometers to kilometers—have mainly been focused on homogenization, replacing the original porous medium by an equivalent spatially structureless one.

The most famous approach to the scale-up problem along these lines is thermodynamically constrained averaging theory (TCAT) [10–14], based on thermodynamically consistent definitions made at the continuum scale based on volume averages of pore-scale thermodynamic quantities, combined with closure relations based on homogenization [15]. All variables in TCAT are defined in terms of pore-scale variables. However, these result in many variables and complicated assumptions are needed to derive useful results.

Another homogenization-approach based on non-equilibrium thermodynamics uses Euler homogeneity to define the up-scaled pressure. From this, Kjelstrup et al. derived constitutive equations for the flow while keeping the number of variables down [16–18].

There is also an ongoing effort in constructing a scaled-up theory based on geometrical properties of pore space by using the Hadwiger theorem [19–21]. The theorem implies that we can express the properties of the spatial geometry of the three-dimensional porous medium as a linear combination of four Minkowski functionals: volume, surface area, mean curvature, and Gaussian curvature. These are the only four numbers required to characterize the geometric state of the porous medium. The connectivity of the fluids is described by the Euler characteristic, which by the Gauss–Bonnet theorem can be computed from the total curvature.

A different class of theories is based on detailed and specific assumptions concerning the physics involved. An example is local porosity theory [22–27]. Another example is the decomposition in prototype flow (DeProf) theory which is a fluid mechanical model combined with non-equilibrium statistical mechanics based on a classification scheme of fluid configurations at the pore level [28–30]. A third example is that of Xu and Louge [31] who introduced a simple model based on Ising-like statistical mechanics to mimic the motion of the immiscible fluids at the pore scale, and then scale up by calculating the mean field behavior of the model. In this way, they obtain the wetting fluid retention curve.

The approach we take in this paper is built on the approach to the upscaling problem found in [32–36]. The underlying idea is to map the flow of immiscible fluids in porous media, a dissipative and therefore out-of-equilibrium system, onto a system which is in equilibrium. The Jaynes principle of maximum entropy [37] may then be invoked and the scale-up problem is transformed into that of calculating a partition function [35].

In order to sketch this approach, we need to establish the concept of steady-state flow. In the field of porous media, “steady-state flow” has two meanings. The traditional one is to define it as flow where all fluid interfaces remain static [38]. The other one, which we adopt here, is to state that it is flow where the macroscopic variables remain fixed or fluctuate around well-defined averages. The fluid interfaces will move and on larger scales than the pore scale where one will see fluid clusters breaking up and merging. If the porous medium is statistically homogeneous, we will see the same local statistics describing the fluids everywhere in the porous medium [36].

Imagine a porous plug as shown in Figure 1. There is a mixture of two immiscible fluids flowing through it in the direction of the cylinder axis under steady-state conditions.

A central concept in the following is the representative elementary area (REA) [34, 35, 39]. We pick a point in the porous plug. In a neighborhood of this point, there will a set of streamlines associated with the velocity field v ⃗ . The overall flow direction is then defined by the tangent vectors to the streamlines. We place an imaginary disk of area A at the point, with orientation such that A is orthogonal to the overall flow direction. We illustrate this as the lower disk in Figure 1. We assume that the disk is small enough for the porous medium to be homogeneous over the size of the plane with respect to porosity and permeability. This disk is an REA.

The REA will contain different fluid configurations. By “fluid configuration”, we mean the spatial distributions of the two fluids in the disk and their scalar velocity fields. We also show a second disk in Figure 1, which represents another REA placed further along the average flow. Since the flow is incompressible, we can imagine the z-axis corresponding to a pseudo-time axis. We may state that the fluid configuration in the lower disk of Figure 1 evolves into the fluid configuration in the upper disk under pseudo-time-translation. We now associate entropy to the fluid configurations in the sense of Shannon [40]. Even though the system is producing molecular entropy through dissipation, it is not producing Shannon entropy along the z-axis. This allows us to use the Jaynes maximum entropy principle, and a statistical mechanics based on Shannon entropy ensues [35].

This statistical mechanics scales up the pore-scale physics, represented by the fluid configurations to the Darcy scale, which then is represented by a thermodynamics-like formalism involving averaged velocities.

Rudiments of this thermodynamics-like formalism was first studied by Hansen et al. [32], who used extensivity to derive a set of equations that relate the seepage velocity of each of the two immiscible fluids flowing under steady state conditions through a representative volume element (REV). They introduced a co-moving velocity together with the average seepage velocity v that contained the necessary information to determine the seepage velocities of the more wetting and the less wetting fluids, v _w and v _n respectively. Stated in a different way, these equations made it possible to make the transformation v , v m ⇆ v w , v n . (1)

It is the aim of this paper to present a geometrical interpretation of the transformation in Eq. 1. We wish to provide an intuitive understanding of precisely what the transformation means and to determine the role of the co-moving velocity. Our aim is not to develop the theory presented in [32] further by including new results, but rather consolidate and amend the existing theory with a deeper understanding.

Our geometrical interpretation will rest upon the definition of a space spanned by two central variables: the wetting and non-wetting transversal pore areas. These areas are defined as the area of REAs covered by the wetting and non-wetting fluids, respectively.

Instead of using the two extensive areas directly, one often works with coordinates where one is the (wetting) saturation. If the pore area of the REA is fixed, it is often convenient to let the other variable be this area. A third way—which is new—is to use polar coordinates. This, as we shall see, simplify the theory considerably.

In Section 2, we describe the REA and the relevant variables. We define intensive and extensive variables, tying them to how they scale under scaling of the size of the REV. We also review the central results of Hansen et al. [32] here: the introduction of the auxiliary thermodynamic velocities and the co-moving velocity.

The central concept in our approach to the immiscible two-fluid flow problem is that of the pore areas. In Section 3 we introduce area space, a space that simplifies the analysis considerably. We showcase different coordinate systems on this space, focusing on polar coordinates.

Section 4 expresses the central results of Hansen et al. [32] in polar coordinates. Expressed in this coordinate system, the geometrical structure implied by the extensivity versus intensivity conditions imposes restrictions that simplify the equations. The section goes on to express the total pore velocity, seepage velocities, thermodynamic velocities, and co-moving velocity in terms of polar coordinates. This provides a clear geometrical interpretation of the co-moving velocity.

Section 5 focuses on relations between the different velocities that can be expressed without referring to any coordinate system. We do this by invoking differential forms and exterior algebra [41, 42]. Differential forms constitute the generalization of the infinitesimal line, area, or volume elements used in integrals and exterior calculus is the algebra that makes it possible to handle them. Flanders predicted in the early sixties [41] that within short they would become important tools in engineering. This did not happen, and they have remained primarily within mathematics and theoretical physics. Such objects are also especially prevalent in thermodynamics, and here, we will show that the computational rules of differential forms provide a simple way of expressing the relations between the velocities in the system at hand.

Section 6 focuses the discussion on the co-moving velocity by constructing a coordinate-independent expression that defines it.

We present in Section 7 two concrete systems. 1) A porous medium obeying the Brooks–Corey relative permeability and 2) a capillary fiber bundle model, which we analyze using the methods developed in this paper.

In Section 8, we summarize our main results which may be stated through three Eqs 58, 69, 98. These equations formulate in a coordinate-free way the three relations that exist between the average seepage velocity and the co-moving velocity.

We will now elaborate on the definition of the REA presented in Section 1.

We define the transversal pore area A _p as previously defined, which defines porosity ϕ = A p A . (2)

The pore volume contains two immiscible fluids; the more wetting fluid (to be referred to as the wetting fluid) or the less wetting fluid (to be referred to as the non-wetting fluid). The transversal pore area A _p may, therefore, be split into the area of the REA disk cutting through the wetting fluid A _w or cutting through the non-wetting fluid A _n so that we have A p = A w + A n . (3)

We also define the saturations S w = A w A p , (4) S n = A n A p , (5) so that S w + S n = 1 . (6)

We note that the transversal pore areas A _p , A _w , and A _n are extensive in the REA A, that is A _p → λA _p , A _w → λA _w , and A _n → λA _n when A → λA. The porosity ϕ and the two saturations S _w and S _n are intensive in the REA A: ϕ → λ ⁰ ϕ, S _w → λ ⁰ S _w and S _n → λ ⁰ S _n when A → λA.

There is a time averaged volumetric flow rate Q through the REA. The volumetric flow rate consists of two components, Q _w and Q _n , which are the volumetric flow rates of the wetting and the non-wetting fluids. We have Q = Q w + Q n . (7)

We define the total, wetting, and non-wetting seepage velocities, respectively, as v = Q A p , (8) v w = Q w A w , (9) v n = Q n A n . (10)

Using Eqs 7–10, we find v = Q A p = A w A p Q w A w + A n A p Q n A n = S w v w + S n v n . (11)

The volumetric flow rates Q, Q _w , and Q _n are extensive in the REA A, and the velocities v, v _w , and v _n are intensive in the area A.

The transversal pore areas A _w ≥ 0 and A _n ≥ 0 parametrize the first quadrant of R 2 . This serves as the “area space” mentioned earlier, which we will continue to denote as such. It is important to realize that this space is not physical space. A value of the transversal pore area A _p corresponds to a point A w , A n in this space. If we consider the entire plane R 2 as a whole and simply restrict our attention to the first quadrant, we may treat the area space as a vector space. The point A w , A n may then equivalently be described by a vector, A ⃗ = A w e ⃗ w + A n e ⃗ n , (30) where e ⃗ w and e ⃗ n form an orthonormal basis set, shown in the leftmost figure in Figure 2. Note that in this picture, the bases are shown as attached to the point A w , A n . We reiterate that both A _w and A _n are extensive.

Every point in the transversal area space corresponds to a given saturation S _w and a transversal pore area A _p and we may view the map (A _w , A _n ) → (A _p , S _w ) given by Eqs 3, 4 as a coordinate transformation. This is a more natural coordinate system to work with since S _w is an intensive variable and A _p is in practice kept constant. We name this system the saturation coordinate system. We show the normalized basis vector set in the middle figure in Figure 2. We calculate u ⃗ p = ∂ A ⃗ ∂ A p S w = S w e ⃗ w + 1 − S w e ⃗ n , (31) u ⃗ s = ∂ A ⃗ ∂ S w A p = A p e ⃗ w − A p e ⃗ n . (32)

We normalize the two vectors u ⃗ p and u ⃗ s to find e ⃗ p = S w e ⃗ w + 1 − S w e ⃗ n 1 − 2 S w 1 − S w 1 / 2 , (33) e ⃗ s = e ⃗ w − e ⃗ n 2 . (34)

This basis set is not orthonormal, as illustrated in Figure 2. By expressing A ⃗ in this coordinate system, we find A ⃗ = A p 1 − 2 S w 1 − S w 1 / 2 e ⃗ p . (35)

As it will become apparent, the most convenient coordinate system to work with from a theoretical point of view is polar coordinates A r , ϕ , given by A r = A w 2 + A n 2 , (36) ϕ = arctan A n A w , (37) or vice versa. A w = A r ⁡ cos ⁡ ϕ , (38) A n = A r ⁡ sin ⁡ ϕ . (39)

We note that ϕ is intensive and A _r is extensive. We show the polar coordinate system in the lower figure in Figure 2.

The basis vector set ( e ⃗ r , e ⃗ ϕ ) is orthonormal, where e ⃗ r = cos ⁡ ϕ e ⃗ w + sin ⁡ ϕ e ⃗ n = e ⃗ p , (40) e ⃗ ϕ = − sin ⁡ ϕ e ⃗ w + cos ⁡ ϕ e ⃗ n . (41)

As for the saturation coordinate system, there is one intensive variable, ϕ, and one extensive variable, A _r . However, in contrast to the saturation coordinate system, both variables are varied in practical situations.

We have that A ⃗ = A r e ⃗ r . (42)

We see that this is consistent with Eq. 35 since A r = A p 1 − 2 S w 1 − S w 1 / 2 (43) and e ⃗ p = e ⃗ r .

We now focus on Eq. 12 which we repeat here. Q λ A w , λ A n = λ Q A w , A n .

We may interpret this equation geometrically. If we follow the value of Q along a ray passing through the origin of the two-dimensional space spanned by (A _w , A _n ) keeping the ratio A _n /A _w constant, it grows linearly with the distance from the Q-axis.

In polar coordinates, this means Q A r , ϕ = v ̂ r A r , (44) where v ̂ r = ∂ Q ∂ A r ϕ = v ̂ r ϕ . (45)

The important point here is that v ̂ r is not a function of A _r .

We may derive Eq. 44 from Eq. 12 as follows. First, we write A _w and A _n in terms of A _r and ϕ in Eq. 12 so that 1 λ Q λ A r ⁡ cos ⁡ ϕ , λ A r ⁡ sin ⁡ ϕ = Q A r ⁡ cos ⁡ ϕ , A r ⁡ sin ⁡ ϕ = Q A r , ϕ . (46)

Next, we set λ = 1/A _r to find A r Q cos ⁡ ϕ , sin ⁡ ϕ = A r v ̂ r ϕ = Q A r , ϕ , (47) where v ̂ r ϕ ≡ Q cos ⁡ ϕ , sin ⁡ ϕ , (48) and we are done. We illustrate in Figure 3 the geometrical meaning of Eq. 12, that is, Euler homogeneity. The graph of the volumetric flow rate Q takes the shape of a “crumpled” cone in the space spanned by (Q, A _w , A _n ).

We now express the different fluid velocities, namely the pore velocity, seepage velocities, thermodynamic velocities, and the co-moving velocity, in terms of polar coordinates.

Since the difference between the thermodynamic- and seepage velocities in the previous section was shown to sit in the tangent vector components, we can benefit from examining the problem from a position where these components are easier to work with. In this section, we therefore introduce differential forms and exterior algebra [41, 42] to make this treatment more economic. In addition, this makes it easier to formulate the equations met so far in a way that does not depend on the coordinate system used.

Eq. 44, Q = v ̂ r A r , may be differentiated to give d Q = v ̂ r d A r + d v ̂ r A r = v ̂ r d A r + v ̂ r ′ A r d ϕ , (70) where v ̂ r ′ = d v ̂ r / d ϕ . Here, dQ, dA _r , and dϕ are one-forms. Furthermore, we have that d Q = ∂ Q ∂ A r ϕ d A r + ∂ Q ∂ ϕ A r d ϕ = v ̂ r d A r + v ̂ ϕ A r d ϕ . (71)

Comparing Eqs 70, 71 gives v ̂ ϕ = v ̂ r ′ . (72)

We do the same using the coordinate system (A _w , A _n ). From Eq. 50, we find d Q = v ̂ w d A w + v ̂ n d A n + A w d v ̂ w + A n d v ̂ n = v ̂ w d A w + v ̂ n d A n . (73)

In the second equality, we assumed that dQ is an exact differential, and in the first equality, we applied the differential to Q = v ̂ w A w + v ̂ n A n . This implies the relation A w d v ̂ w + A n d v ̂ n = 0 . (74)

This equation corresponds to the Gibbs–Duhem equation in thermodynamics, which is a statement of the dependency amongst the intensive variables.

We express v ̂ w and v ̂ n in terms of v ̂ r and v ̂ ϕ by inverting Eqs 61, 62 to find v ̂ w = v ̂ r ⁡ cos ⁡ ϕ − v ̂ ϕ ⁡ sin ⁡ ϕ , , (75) v ̂ n = v ̂ r ⁡ sin ⁡ ϕ + v ̂ ϕ ⁡ cos ⁡ ϕ . (76)

By applying the exterior derivative, we get the relations d v ̂ w = cos ⁡ ϕ d v ̂ r − sin ⁡ ϕ d v ̂ ϕ − v ̂ n d ϕ , (77) d v ̂ n = sin ⁡ ϕ d v ̂ r + cos ⁡ ϕ d v ̂ ϕ + v ̂ w d ϕ . (78)

Combining Eqs 75–78 with Gibbs–Duhem Eq. 74 gives d v ̂ r d ϕ = v ̂ ϕ , (79) which is identical to Eq. 72. Hence, this is the Gibbs–Duhem equation in polar coordinates.

The thermodynamic velocity, Eq. 59, is, therefore, given by v ⃗ ̂ = v ̂ r e ⃗ r + v ̂ r ′ e ⃗ ϕ . (80)

We now have one of the main results of Hansen et al. [32] in a very compact form. v ⃗ = v ̂ r e ⃗ r + v ̂ r ′ + v ̂ m e ⃗ ϕ . (81)

We note that since v ̂ r only depends on ϕ, the same must be true for v ̂ m . Hence, v ̂ m = v ̂ m ( ϕ ) . Here, we see that the relation between the thermodynamic- and seepage velocities are not captured in the extensive structure tied to the radial coordinate, but instead in the intensive quantities which only has a dependency on ϕ.

We now differentiate Eq. 56, Q = v ⃗ ⋅ A ⃗ written in the saturation coordinate system, to find d Q = v d A p + d v d S w A p d S w . (82) We note that d A p = d A w + d A n , (83) from differentiating Eq. 3. Likewise, we note that d S w = ∂ S w ∂ A w A n d A w + ∂ S w ∂ A n A w d A n = 1 A p 2 A n d A w − A w d A n , (84) using Eq. 4. We combine Eqs 82–84 and find d Q = v + 1 − S w d v d S w d A w + v − S w d v d S w d A n . (85)

By using Eq. 13, definitions of Eqs 15, 16, and transformations Eqs 19, 20, we may write d Q = v ̂ w d A w + v ̂ n d A n = v w − 1 − S w v m d A w + v n + S w v m d A n . (86)

Equating Eq. 85 and Eq. 86 gives Eqs 24, 25.

Let us now rewrite the differential dQ as follows: d Q = d A p v = d A p S w v ̂ w + 1 − S w v ̂ n = A p v ̂ w − v ̂ n d S w + v d A p , (87) where we have used Eq. 74. We now combine this equation with Eq. 82 to get d v d S w = v ̂ w − v ̂ n . (88)

By using Eqs 19, 20, we may rewrite this equation in terms of the seepage velocities, resulting in Eq. 26.

We differentiate Q a second time, which is zero according to the Poincaré lemma [41, 42]. We see that this is indeed so by using Eq. 70, d 2 Q = d v ̂ r d A r + v ̂ r ′ A r d ϕ = v ̂ r ′ d ϕ ∧ d A r + v ̂ r ′ d A r ∧ d ϕ = v ̂ r ′ − v ̂ r ′ d ϕ ∧ d A r = 0 , (89) where the wedge ∧ signifies the antisymmetric exterior product.

The same equation expressed in the coordinate system (A _w , A _n ) gives d 2 Q = d v ̂ w d A w + v ̂ n d A n = ∂ v ̂ n ∂ A w A n − ∂ v ̂ w ∂ A n A n d A w ∧ d A n = 0 . (90)

We rewrite the coefficient in terms of the saturation coordinate system, ∂ v ̂ n ∂ A w A n − ∂ v ̂ w ∂ A n A n = d v ̂ n d S w ∂ S w ∂ A w A n − d v ̂ w d S w ∂ S w ∂ A n A w = S w d v ̂ w d S w + 1 − S w d v ̂ n d S w = 0 , (91) which is nothing but the Gibbs–Duhem equation yet again. Writing this equation in terms of the seepage velocities gives us Eq. 28.

We see that the different equations relating the velocities can all be found from the differential geometric structure of the space spanned by the areas, with addition to Euler homogeneity, see Figure 3. All relations are either a consequence of writing dQ using different coordinate systems or a consequence of the Poincaré lemma, expressed as d ² Q = 0.

The co-moving velocity appears in several equations in the previous sections. Common to all of them is that they all are written out in terms of a given coordinate system. We may write the representation of the co-moving velocity in the saturation coordinate system, Eq. 28, as A p v m = A ⃗ ⋅ d v ⃗ d S w . (92)

We may write Eq. 65 in polar coordinates as A r v ̂ m = − A ⃗ ⋅ d v ⃗ d ϕ . (93)

Expressed in terms of the seepage velocities, we have that d Q = v w d A w + v n d A n + A w d v w + A n d v n . (94)

In order to obtain an expression for the co-moving velocity which is independent of the coordinate system, we construct A w d v w + A n d v n = v m A n d A w − A w d A n A w + A n = v m A p d S w = − v ̂ m A r d ϕ , (95) where the second line only reflects Eqs 92, 93.

We write dQ as d Q = d [ A ⃗ ⋅ v ⃗ ̂ ] = v ̂ i d A i , (96) using Eq. 74. Here, Einstein summation convention has been applied, with the index i running over w , n in the basis A w , A n . In terms of the seepage velocities, this becomes d Q = d [ A ⃗ ⋅ v ⃗ ] = v i d A i + A i d v i . (97)

Combining Eqs 96, 97, 66 gives us v ̂ m i d A i + A i d v i = 0 , (98) where v ̂ m i denotes the components of v ⃗ ̂ m . This holds in any coordinate system, since it is just obtained by differentiating a function, namely A ⃗ ⋅ v ⃗ . We have, therefore, obtained an expression for the co-moving velocity which is independent of the chosen coordinate system.

We can introduce vector-valued differential forms [42] to make a connection between the interpretation in terms of (tangent) vectors used up until now and the language of differential forms. Let e ⃗ i be an arbitrary basis for vectors on the area space as we have discussed until now, be it area-, saturation-, or polar coordinates. We now take A ⃗ as an example. We write A = β i e ⃗ i , (99) where we have applied the Einstein summation convention. A is now regarded as a vector-valued differential form (thus seeing e ⃗ i as components written out in the β ⁱ basis, i.e., a reversal of viewpoint), the β ⁱ are 0-forms, just functions, and the index i runs over the coordinate indices. We can still treat A as a vector, but the difference in viewpoints is apparent when we apply the exterior derivative d to A, which gives d A = ∑ i e ⃗ i ⊗ d β i = d β i e ⃗ i , (100) where we have suppressed the tensor product of the forms and basis vectors since we are working over a vector space, as is standard notation in for example [42]. Strictly speaking, the exterior derivative d should be replaced by the exterior covariant derivative in this setting, a generalization of d which is defined on both tangent vectors and forms (there is currently no term in Eq. 100 that reflect the changes in the basis e ⃗ i ). The terms A ⁱ dv _i are then possible to interpret in terms of connection forms. However, this is outside the scope of the current discussion, and we leave this topic for future work.

If we write out Eq. 100 in the coordinate system A w , A n , we get d A = d A w e ⃗ w + d A n e ⃗ n . (101)

With this formalism in mind, we now return to polar coordinates. Eq. 98 may be written in these coordinates by using Eqs 40, 41 to define the vector valued 0-forms, ω r = cos ⁡ ϕ e ⃗ w + sin ⁡ ϕ e ⃗ n , , (102) ω ϕ = − sin ⁡ ϕ e ⃗ w + cos ⁡ ϕ e ⃗ n , (103) which are equivalent to the corresponding basis vectors in Eqs 40, 41, but interpreted as 0-forms. We then compute d ω r = e ⃗ ϕ d ϕ , , (104) d ω ϕ = − ω r d ϕ (105) so that d A = d A r ω r = d A r ω r + A r d ω r = d A r e ⃗ r + A r e ⃗ ϕ d ϕ , (106) where in the second equality, we have used ω _r which is just another way of writing the basis vector e ⃗ r that emphasizes its role as a form.

We can similarly define vector valued 1-forms d v ̂ and dv from Eqs 59, 49 respectively. We then straightforwardly obtain d v ̂ = d v ̂ w e ⃗ w + d v ̂ n e ⃗ n , , (107) d v = d v w e ⃗ w + d v n e ⃗ n . (108)

We express dv in polar coordinates as d v = v r ′ e ⃗ r d ϕ + v r e ⃗ ϕ d ϕ + v ϕ ′ e ⃗ ϕ d ϕ − v ϕ e ⃗ r d ϕ = v r ′ − v ϕ e ⃗ r + v r + v ϕ ′ e ⃗ ϕ d ϕ , (109) giving A r v ̂ m d ϕ = A r v ϕ − v r ′ d ϕ , (110) or v ̂ m = v ϕ − v r ′ . (111)

By combining Eqs 54, 68, 72, we obtain the same equation. Eq. 98 written in the saturation coordinate system yields Eq. 28.

We give in this section a couple of examples of practical use of the ideas that have been presented. We consider first the Brooks–Corey relative permeability model [43] and then a capillary fiber bundle model [44, 45].