Geometrically derived efficiency of slow immiscible displacement in porous media

The efficiency of a displacement is the fraction of applied work over the change in free energy. This displacement efficiency is essential for linking wettability to applied work during displacement processes. We quantify the efficiency of slow immiscible displacements in porous media from pore space geometry. For this end, we introduce pore-scale definitions for thermodynamically reversible (ison) and irreverisble (rheon) processes. We argue that the efficiency of slow primary displacement is described by the geometry of the pore space for porous media with a sufficient number of pore bodies. This article introduces how to calculate such geometry-based efficiency locally, and integrating this local efficiency over the pore space yields an aggregate efficiency for the primary displacement in the porous medium. Further, we show how the geometrical characterization of the displacement efficiency links the efficiency to the constriction factor from transport processes governed by the Laplace equation. This enables estimation of displacement efficiency from traditional and widely available measurements for porous media. We present a thermodynamically based wettability calculation based on the local efficiency and a method to approximate this thermodynamically based wettability from traditional experiments.


I. INTRODUCTION
Wettability is the relative adherence of different fluid phases to the substrate. For a system of two fluid phases, the phase with the largest affinity to the solid is called the wetting phase, while the other is called the nonwetting phase. In the first modern account of capillary action, going back more than 200 years, Leslie noted that heat is liberated when water or oil is added to paper or linen [1]. Paper and linen are porous media, and air has a larger interfacial tension with the solid surface inside these porous media than both water and oil. The interfacial tension can be defined in terms of surface energy per area, thus paper and linen filled with air have a higher energy than when they are filled with water or oil. This energy difference leads to liberation of heat when water or oil is imbibing. Two decades after Leslie the link between interfacial tensions and heat of wetting was further quantified when different temperature increases was reported for a range of fluid and solid combinations [2].
The wettability between two fluid phases can be regarded as a force balance between the interfacial tensions, a perspective taken by Young a couple of years after the seminal work of Leslie [3]. For a flat solid surface, the force balance in direction of the surface includes the contact angle, i.e., the angle between the fluid-fluid surface and the solid surface at the three-phase contact line. This contact angle is considered a measure of wettability. In contrast to the force-balance perspective of Young, wettability can be viewed as an energy balance of free energies of interfaces, a slightly newer view tracing back to Gauss [4].
The energy-balance equation of Gauss can be extended by including the free energy of three-phase contact lines [5], as suggested in a footnote by Gibbs [[6], p. 455]. If one assumes that the Hadwiger theorem holds for free energy in porous media, then the free energy is described as a linear combination of intrinsic volumes of the fluid and solid phases [7,8]. This gives a free-energy description by fully 27 variables, fortunately geometrical dependencies reduce the number to seven independent geometrical descriptors consisting of the interfacial areas, three-phase contact lines, interface curvatures and the Euler characteristic [9]. However, the interfacial energies are often considered dominant in porous media with μm-size pores. In the following we will therefore equate the free energy with the interfacial energies, assuming that other contributions to the free energy are of minor importance.
When equating the free energy of a system consisting of two fluid-phases inside a porous medium with the interfacial energies, the work applied to the system during a reversible displacement process can be equated with changes in the interfacial energies [10,11]. When one incompressible phase is entering the porous medium, displacing another incompressible phase leaving the medium, the pressure difference between the entering and produced fluids times the volume change of the fluids yields the work applied to the system. When equating this applied work to the changes in interfacial energies, several authors have estimated the wettability from measurements of applied work and changes in interfaces [9,12,13].
The applicability of using this thermodynamic relationship to infer wettability in porous media has been hampered by the assumption of a reversible displacement. The efficiency of the displacement, i.e., the fraction of change in interfacial energy to applied work, is dependent on flow rate and pore structure [10,14]. Except for simplified pore structures and slow flow, the displacement is nonreversible, yielding an efficiency lower than one. Nonreversible pore-scale processes induce pressure signals that have been interpreted to infer pore structure and wettability information [15,16]. However, traditional empirical wettability measurements, such as the Amott and USBM tests, do not adjust for pore structure effects [17]. For these methods the pore structure is then a determining factor for the measured wettability. Recent developments in imaging enable description of the pore structure at a resolution which makes possible reliable computation of transport properties and displacement processes [18][19][20]. This has opened up new possibilities for describing transport properties of porous media by geometrical descriptors [21]. In this article we will investigate the correspondence between pore structure and efficiency, thereby untangling pore structure from wettability.
Displacement and transport in porous media is of importance for a wide range of industrial processes, including fluid flow in partially saturated soils and displacement of hydrocarbons by injected water [22]. Wettability has a first-order effect on such multiphase flow. In previous work [10,23] the efficiency of displacement processes has been estimated based on images of interface areas during displacement. In the present work we show that the efficiency of slow immiscible displacement can be obtained from the static pore structure through a pore geometry description. Our derivations will yield a local efficiency, which is applied to a set of three-dimensional porescale images and models to predict efficiency as a function of saturation. Further, the derived efficiency is applied to an imaged displacement experiment for estimation of wettability. We also investigate cross-correlations between our developed pore-structure-based displacement efficiencies and other porestructure descriptors, which open up the possibility to estimate displacement efficiency based on more readily available measurements.

II. THERMODYNAMICALLY REVERSIBLE AND IRREVERSIBLE DISPLACEMENT PROCESSES
In this section we will introduce pore-scale definitions for thermodynamically reversible and irreverisble processes. These definitions will be used in the remainder of the article. Assume a porous medium V of length s in direction of an applied pressure difference p. The porous medium V consists of matrix and pore space ⊂ V , where the pore space is filled with two immicible incompressible fluids. For our considerations we assume that the length scales are small enough to disregard gravity, however, all derivations in this article can be extended to include gravity. We assume that the interfaces between the two fluids are infinitesimally thin, and we allow for discontinuity of the pressure field at the interface. The fluids might have different viscosity μ i (and density ρ i , only of importance for our calculations when gravity is included), where the index i = 1, 2 indicates the two different fluids. We assume a slow flow, where the interstitial velocity field u is governed by the Stokes equation with a no-slip condition at the fluid-solid interface.
In our system an external total flow rate Q is imposed on the system. This flow rate induces the pressure difference p = p out − p in over the porous medium, where we assume that the inlet pressure p in and outlet pressure p out are constant over their respective end-face for any given time. When the total flow rate is Q and the pressure difference is p, then the rate of applied energy is − pQ (note that p is negative). The potential driving the interstitial fluid flow is the interstitial pressure gradient ∇ p, and the rate of work done by this pressure gradient is given by −∇ p · u, where u is the interstitial velocity. In where n is the outward pointing unit normal field of the fluidfluid boundaries δ f . Note that these fluid boundaries δ f are on both sides of the fluid-fluid interfaces, with one unit normal pointing out of fluid 1 and the corresponding pressure inside fluid 1, while the inverted unit normal is pointing out of fluid 2 with pressure in fluid 2. Thus, each fluid-fluid interface gives two fluid boundaries. The left side of Eq. (1) is the applied power, while the first term on the right side equals the viscous dissipation inside the pore space. The last term yields the rate of change in (Helmholtz) free energy in the system, i.e., the rate of change in free energy from the interfacial tensions (and possibly other energies, e.g., energy from the contact lines). The free energy associated with the placement of an individual fluid-fluid interface can either increase or decrease. Let us start by considering the case where the free energy increases for all internal interfaces. Then the interstitial fluid flow is driven by the externally applied pressure drop only. Consider an infinitesimal external flow rate Q → 0 (slow flow). Due to linearity in the Stokes equation, ∇ p → 0 and u → 0 at the same rate as Q → 0, while the pressure drop over the fluid-fluid interfaces can result in a p 0 when Q → 0. Then for the nontrivial case (i.e., p 0 when Q → 0) we see that ∇ p · udV goes to zero faster than pQ, thus we get pQ ∇ p · udV when Q → 0. Rearranging Eq. (1) then gives Thus, in this slow-flow limiting case viscous dissipation vanishes, and the externally applied energy is expended on increasing the free energy related to the placement of the fluid-fluid interfaces. For higher flow rates Q, some energy will be expended on viscous flow in addition to the increase in free energy of the surfaces. This will give a nonzero second term in Fig. 2, consisting of nonreversible viscous dissipation of energy. Thus, the assumption of a reversible displacement only holds for slow flow. In the following we will assume that the flow is slow enough, that the wettability is constant and without sub-resolution heterogeneity and that the solid surfaces are smooth, such that the similarity in Eq. further discussion on sub-resolution heterogeneity and surface roughness.
We define a displacement where as an ison, loosely following the terminology in Ref. [10]. In other words, when the similarity in Eq. (2) holds, we have an ison displacement. An ison displacement is a reversible process. During an ison displacement one can use the thermodynamic relationship between applied work and changes in interfaces to infer the wettability.
We will now turn to irreversible processes. For the reversible case the change in free energy was positive for all fluid-fluid interfaces. Contrary to the reversible case, assume that there exists at least one fluid-fluid interface releasing energy. Let {δ i f | i ∈ S} be the set of individual fluid-fluid interfaces. Thus, assume there exists an i ∈ S such that the integral of p u · n over the two sides of the interface δ i f is negative; δ i f p u · ndS < 0. This release of energy will induce fluid flow, thus u 0. As we consider a slow displacement, Q 0. Thus, the induced flow leads to pQ ∇ p · udV , which implies that With infinitesimal external flow, Q 0, we can assume that no external work is applied to the system during the time span of an internal redistribution of fluids. Thus, the excess free energy of the system must decrease during this internal Thus, for a displacement where movement of one fluid-fluid interface leads to release of energy, then the movement of all the fluid interfaces will lead to a decrease in the total free energy of all the interfaces. We will define a rheon displacement as a displacement where again loosely following the terminology in Ref. [10]. A rheon displacement is thus a thermodynamically irreversible internal movement of interfaces at a constant saturation. Haines jumps are examples of rheon displacements [14]. The energy expended during a rheon displacement is extracted from the interfacial energies, thus this dissipated energy must be subtracted from the applied work when estimating the wettability from the correspondence between applied work and interfacial energies. This will be applied in Sec. VIII.

III. TWO-PHASE DISPLACEMENT IN SIMILAR CIRCULAR PORE TUBES
In this section we will consider two-phase primary displacement in bundles of pore tubes to illustrate how the amount of fluid-fluid interfaces affect the efficiency of the displacement process. By primary displacement, we mean the displacement process which starts with the porous medium filled by a single fluid and then displaced by another immiscible fluid. In all our examples only one fluid, i.e., the displaced fluid, will leave the medium at one side, while the injected fluid will be injected at the opposite side.
We will now consider a set of tubes with varying crosssectional area. As an example case, these tubes are 200 μm long with a sinusoidal shape, where the largest cross-sectional radius is 20 μm while the smallest radius is 10 μm, as depicted in Fig. 1. We assume that the tubes are initially filled with fluid 1, while connected to fluid 1 at the right side and fluid 2 at the left side. Fluid 1 is wetting the solid, thus the solid-fluid interfacial tension between fluid 2 and the surface, σ 2s , is assumed to be larger than between fluid 1 and the surface; σ 2s > σ 1s . Thus, there is a higher excess free energy of the system when fluid 2 is in contact with the surface than when fluid 1 is in contact with the surface.
As the flow is infinitesimally slow the pressure in each connected fluid phase can be considered constant during the reversible isons. Further, the fluid-fluid interfaces can be considered at equilibrium, so the pressure difference between the two fluids is assumed to be given by the Young-Laplace equation, where σ is the fluid-fluid interfacial tension, while κ = 1/r 1 + 1/r 2 is the mean curvature, with r 1 and r 2 being the principal radii of curvature. As the cross-section of our sinusoidal tubes are assumed to be circular, the principal radii of curvature will be equal, say r 1 = r 2 = r, thus the mean curvature is κ = 2/r. Initially there is no pressure difference between the two fluids, corresponding to a flat fluid-fluid interface at the left inlet (κ = 0). We will inject fluid 2 at an infinitesimal but constant rate Q. This will give the surface a spherical shape, where the radius of curvature is decreasing with increasing pressure. This pressure increase is seen at the right side of the pressure plots in Fig. 2.
Eventually, the contact line between the fluid-fluid interface and the solid will start moving into the sinusoidal tube. As we are moving the fluids at an infinitesimal rate, and therefore are assuming equilibrium conditions, we assume that the contact angle between the fluid-fluid interface and the solid surface is given by Young's equation: where σ is is the interfacial tension between the two fluids i = 1, 2 and the solid s, while σ d = σ 2s − σ 1s is the difference in fluid-solid tensions. Note that if σ 2s − σ 1s = σ d > σ , then it is more energetically beneficial to leave a thin layer of fluid 1 to coat the surface after fluid 2 has filled the bulk of the tube, and we can assume that the contact angle is zero. One could consider the fraction ω i = σ d /σ as a wetting index in its own right, with this wetting index being indicative of changes in wetting properties also in the extreme cases where ω i > 1.
When the three-phase contact line is at the smallest constriction, as indicated by the full-drawn red line in Fig. 1, then the fluid configuration is unstable. An infinitesimally larger cross-sectional radius for one tube (or any other disturbance to the system) will lead to the fluid-fluid interface moving forward in this tube, while the interfaces in the remaining tubes will move backward. This is a Haines jump, and it is a rheon displacement according to the definition given by Eq. (6). The end-state of the Haines jump has a lower free energy than the initial state. This reduction in excess free energy will partly be expended on viscous fluid flow between the tubes where the interface is moving backward and the tube where the interface is moving forward.
The resulting fluid configuration after the first Haines jump is illustrated by the dashed lines in Fig. 1. We observe that the resulting interfaces have a smaller curvature in the threetube case than in the ten-tube case, which is indicative of the larger pressure drop during the Haines jump in the three-tube case. This is reflected in the larger pressure drops at constant saturation in the left plot versus the right plot in Fig. 2, i.e., the vertical parts of the pressure curve is longer in the left plot than in the right plot. Figure 2 also contains the change in free energy F = F − F i from the initial free energy F i , calculated as where A i is the initial fluid-fluid interfacial area, i.e., the sum of the leftmost cross-sectional areas of the tubes. Note that we do not associate any energy to the three-phase contact line nor the curvature, in contrast to, e.g., Ref. [9]. In Fig. 2 we have also plotted the applied work calculated as 033113-4 where we integrate over the injected volume V 2 of fluid 2. When all the fluid-fluid interfaces move from the inlet towards the outlet, we have since p = p 1 − p 2 and the last integral of the fluid-fluid interfaces gives Q. Thus, according to our definition given by Eq. (3), this is an ison displacement. During the ison displacement the external work equals the change in free energy, dW = dF . Observe that the applied work and free energy curves in Fig. 2 are identical until the first rheon displacement, indicating that this part of the displacement process is a reversible ison displacement. The two curves have a constant separation for all subsequent ison displacements. The applied work stays constant during the Haines jump, as the change in volume of fluid 2 is zero. During a single rheon displacement, the distance traversed by the backward moving fluid interfaces are shorter for the ten-tube case than for the three-tube case. Thus, the pressure difference between the progressing and regressing interfaces are larger for the tentube case than the three-tube case. The loss of free energy will therefore be larger for the ten-tube case than the three-tube case. This is visible when comparing the two plots in Fig. 2.

IV. PRIMARY DISPLACEMENT EFFICIENCY
The efficiency of an immiscible two-phase primary displacement in a porous medium relates how the externally applied work W driving the displacement divides itself into a change in free energy F and energy dissipation (heat production). This can be written as where the displacement efficiency η is expressed as the fraction of the change in free energy F to the externally applied work to the system W [10]. Let us consider different numbers of sinusoidal tubes from those depicted in in Fig. 1. In the case of a single tube, there will be no Haines jumps. As we assume no surface roughness and constant wettability, the displacement process will be fully reversible, i.e., η = 1 (see Appendix A). For an example with two tubes, the movement of the interfaces forward from the first constriction will be reversible, due to the increase in free energy associated with the increase in interfacial area when the interfaces move from the constricted pore throat into the pore bodies.
For three or more tubes the rheons are associated with a reduction in free energy, and part of this energy is expended on the viscous fluid flow during the redistribution of the fluids. As evident from Figs. 1 and 2, a higher number of tubes leads to a larger pressure difference between the set of backward moving interfaces and the single forward moving interface, and thereby a larger drop in excess free energy. Thus, a higher number of tubes will decrease the efficiency of the displacement process. For a high number of tubes the backward moving interfaces will move only slightly, thus we might approximate the backward moving interfaces as stagnant, residing at the smallest constrictions. By the Young-Laplace equation, Eq. (7), combined with Young's equation, Eq. (8), the displacement of fluid 1 by fluid 2 in the larger pore after the first constriction then occurs at a pressure difference of p = p 1 − p 2 = 2σ d /r 2σ d /r t , where r t is the cross-sectional radius of the tube at the pore throat. A reversible ison displacement of the same pore body would have occurred at p = 2σ d /r for the varying cross-sectional area r of the tube. Thus, the displacement efficiency in the setting of a high number of tubes will be For the last equality we have used the earlier assumptions that the fluid-solid interfacial tension difference σ d is constant.
We thus see that for a homogeneous fluid-solid property, this property does not affect the efficiency of the displacement process. In general, when the fluid-solid interfacial tension difference σ d is uncorrelated with the radius r, Eq. (13) holds for sufficiently large samples. The efficiency in the limiting case of a large number of tubes is thus a purely geometrical property. Note that this equation is not restricted to our sinusoidal model, any set of tubes with varying cross-sectional area and where the fluid-solid interfacial tension difference σ d is uncorrelated to the cross-sectional area will be described by the geometrical property given by Eq. (13).
Calculating Eq. (13) for our set of sinusoidal tubes gives a theoretical limiting efficiency of η = 0.7188. The simulated efficiency for different number of tubes are shown in Fig. 3. We observe that the simulated efficiency converges to the theoretical efficiency when the number of tubes increases.
Already around 100 tubes the simulated efficiencies are similar to the theoretical efficiency. This indicates that for a 033113-5 medium where the displacement front covers more than 100 individual pores, one can expect the displacement efficiency η for an event filling a single pore body to be close to the geometrical efficiency given by Eq. (13). While many pores bodies can be filled during a single rheon, natural porous media have fluid-fluid interfaces residing in small pores, in crevasses and in pendular rings, increasing the fluid-fluid area. This large fluid-fluid area sustain a small pressure drop during the rheons.
The transition to smaller pressure drop with increasing system size is supported by experiments and pore-scale simulations. For porous media with a small number of pores the rheons are clearly observable on a plot of the pressure difference curve; see, e.g., Ref. [24]. For slightly larger rock samples, still even as small as in the millimeter range, one needs experimental setups with an extreme resolution to observe the rheons [25]. This indicates that the geometry-based efficiency given by Eq. (13) will be a good approximation for settings such as soil and hydrocarbon reservoirs. Conversely, efficiencies measured or calculated on too small sample sizes are not representative for the displacement efficiencies in pore structures consisting of a high number of pore bodies, e.g., soil and hydrocarbon reservoirs. Now consider a pore structure large enough so that the pressure drop during rheons can be neglected. We continue to assume constant fluid-fluid interfacial tension σ and constant fluid-solid interfacial tensions, where σ d is the difference between the two fluid-solid interfacial tensions as given in Eq. (8). An infinitesimal displacement right before a rheon is assumed to be an reversible ison, thus the pressure difference between the two fluids momentarily before the rheon is given by where the subscript r indicates that this is the pressure difference during the subsequent rheon. Since we have assumed smooth surfaces, the change in fluid-fluid area is considered second order versus the change in fluid-solid area momentarily before a rheon. We thus get the following approximation: Employing Eq. (8), the fraction can be determined by geometrical considerations. For the sinusoidal tubes in Fig. 1 we then have p r /σ d 2/r t , where r t is the cross-sectional area of the pore-throat. Let V r be the volume filled by fluid 2 during a single rheon. Locally, the work done to fill this volume is W r = p r V r , while the local change in free energy during the displacement of volume V r is where A r 12 and A r 2s is the change in the surfaces when fluid 2 fills V r . The efficiency of filling V r is then given by where we have employed p r /σ d 2/r t . In the extreme wetting case, where cos(θ ) = 1 so that σ = σ d , then all the components of F r can be determined geometrically. When A r 12 is small compared to A r 2s there is only a weak dependency between the efficiency and wetting. As an example, for the sinusoidal tubes we have A 12 = 0, thus and the efficiency η r is independent of wetting. We expect A r 12 to be small compared to A r 2s for primary displacement in a range of pore space geometries.
For a volume filled by an ison the efficiency is 1, corresponding to a reversible process. We can thus consider a local efficiency η l , where η l = η r for the rheon volumes V r , while η l = 1 elsewhere. Integrating the local efficiencies will then give the displacement efficiency as The local efficiency can be calculated for any subspace ⊂ , e.g., the subspace spanned by the invading fluid at different saturation steps during an invasion process.

V. GEOMETRICALLY DERIVED PRIMARY DISPLACEMENT EFFICIENCY
We will now turn to complex geometries given by natural porous media and calculate the displacement efficiency based on the local efficiency η l as given in Eq. (20). For simulating the primary drainage process, we will use a pore-morphologybased simulation as introduced in Ref. [26], and later modified in Ref. [27]. For simplicity, we will assume perfect wetting as given by cos(θ ) = 1. The methodology can be extended to other contact angles as in Ref. [28]. There are some wellknown drawbacks with this pore-morphology methodology; the main drawback is that it assumes that the two principal curvatures of the fluid-fluid surface are equal and not infinite, which leads to errors in the Laplace pressure, especially for wetting phase pendular rings between spherical grains. Despite known weaknesses, the methodology was chosen for its efficiency and accuracy for most of the drainage process [27]. The method is based on some basic morphological operations, which can be found in Appendix B.
To calculate the fluid distribution at the different pressures [equivalent to Eq. (B3) in the Appendix], we start by calculating the distance transform d t , giving the Euclidean distance all calculations can be conducted in three dimensions (as they will be in the remainder of this article), the example case used in this section, shown in Figs. 4-6, was restricted to a two-dimensional example for better visualization. The distance transform d t for this simple two-dimensional model is illustrated in Fig. 4.
We then calculate a critical path radius r p (x) for a path from the inlet to point x, that is the minimal distance d t along the path with the largest minimal distance: where S is the set of paths from x to the inlet. Note the similarity with the percolation threshold used in Ref. [29], giving a characteristic length for permeability prediction. We solved for r p with a slightly modified version of Dijkstra's algorithm [30], with Fig. 5 showing the resulting critical path radii for the two-dimensional example case.
We can now look at the balls (or disks in two dimensions) spanned by the value of the distance transform and the critical radii. A point x is assigned to the largest ball covering it as These values can be related to morphological operators as (see The local displacement efficiency, i.e., the relative difference between the irreversible and reversible displacement, is then given by We observe that the local efficiency η l (x) is always smaller or equal to one, as c r p (x) c dt (x). A plot of the local displacement efficiency is shown in Fig. 6.
In Fig. 6 we observe that the local efficiency η l is highest at the inlet side and in the smaller pore tubes. This holds in general; since d t r p for small values of d t , we see that c dt c r p , and thereby η l 1, for small values of d t . Thus, points with small d t values are expected to give high efficiency, hence we expect high efficiency for pore throats, pendular rings, in small crevasses, in the cavities of rough surfaces and other small dead-end pores.

VI. DISPLACEMENT EFFICIENCY FOR DIFFERENT POROUS MEDIA
Using the methodology outlined in Sec. V, we have calculated the local displacement efficiency η l given by Eq. (25) for a set of three-dimensional models and images of natural porous media. The models are Fontainebleau sandstone at different porosity, as used in Ref. [31]. The images are collected from the website of the Imperial College PERM group [32], and consist of both sandstone, sand pack and carbonate samples. Additionally, we have one image of sintered soda lime beads used in Refs. [9,33], which we will return to in the next section.
The calculated local efficiency values for one of the threedimensional Fontainebleau sandstone models are displayed in Fig. 7, with the middle cross-section shown in Fig. 8(a). We observe the same overall effects in this realistic structure as in the simple two-dimensional example in the previous section; close to the inlet and in small pores and crevasses we have a higher efficiency. The differential pressure (capillary pressure) between the two immiscible phases for different saturation values are shown by the colored curve in Fig. 8(c). The colors correspond to the local efficiency for the different saturation values. The black curve corresponds to the displacement efficiency of the drainage process up to the given saturation. The U-shape of the black curve and the colors of the pressure-curve indicate that the efficiency is highest at the beginning (i.e., high saturation values) and at the end (i.e., low saturation values) of the drainage process. As already mentioned, at the end of the drainage process we are filling up the smaller pores and crevasses with high efficiency. At the beginning we are filling pores close to the inlet, an end effect with high efficiency.
The lower efficiency at the inlet, as reflected by the high efficiency values for low saturation values in Fig. 8(c), is an end effect that will affect our effective displacement efficiency for small sample sizes. We will therefore introduce a method correcting for the end effect close to the inlet. For this end, the porous medium can be extended in the direction of applied pressure difference by mirroring. The pressure difference when entering the second copy (and all subsequent copies) of the porous medium will be given by the maximal critical path value at the outlet of the original medium. We will denote the pressure associated to the maximal critical path value at the outlet as the percolation pressure. As indicated by the name, this is the pressure needed for the invading phase to percolate the sample in the direction of the applied pressure. The end effect of high efficiency at the inlet can be removed by assuming that the smallest pressure difference when fluid 2 is entering the medium is given by the percolation pressure. Such a calculation is shown in Fig. 8(b). We observe that there are no end effects in the upper part of the plot, in contrast to Fig. 8(a). The corresponding capillary pressure and efficiency curves are shown in Fig. 8(d). We see the lower efficiency and higher pressure difference for the high saturation values compared to the calculations with end effect shown in Fig. 8(c). For this sample the displacement efficiency is η = 0.65 when using a percolation threshold, compared to a value of 0.67 with the end effect. Thus, this sample is sufficiently large to suppress the influence of the end effect. In general, the calculated displacement efficiency when the starting pressure is given by the maximal percolation path is expected to be more representative for large scale porous media, such as soils and oil reservoirs, as it reduces small size effects. In the remainder of this article all the calculated efficiencies will use the percolation pressure as the minimum injection pressure. However, all samples where large enough to ensure that the end effect had minor impact.
In Fig. 9 we have plotted the displacement efficiency η given by Eq. (20) for the considered models and images. One can observe a trend in the plot, where lower-porosity samples have a lower efficiency. Large pore bodies give large values for c dt , while small pore throats give small c r p values. Thus, from Eq. (25) we see that a large aspect ratio between pore bodies and pore throats leads to low efficiency. The trend seen in Fig. 9 can be linked to larger aspect ratios between pore bodies and pore throats for the lower porosity samples. There is a clear trend for the Fountainebleau samples, which is indicative of a single rock type, created by the same rock forming process [31]. The other samples are more scattered, still there is an overall trend of lower porosity yielding lower efficiency.

VII. CROSS-PROPERTY RELATIONS
As mentioned in the end of the previous section, the efficiency is linked to the aspect ratio between pore bodies and 033113-8  Fig. 7. We note the end effect as low pressures and high efficiencies for high saturation values. Finally, (d) capillary pressure curve when using the percolation threshold for the slow drainage, as visualized in (b). In contrast with (c) we observe lower efficiencies and higher pressures for the high saturation values. pore throats. It is therefore assumed that efficiency would have a stronger correlation with measures for this aspect ratio than with porosity. Slow flow through a cylindrical tube, as given by the solution to the Stokes equations in Eq. (A1), is proportional to the cross-sectional area squared, thus permeability does not have a linear correspondence to the body-throat aspect ratio. A pore-structure description of permeability k as was derived in Ref. [31], where τ s is a streamline-based tortuosity, L h is a characteristic length for the porous medium, while C s is a constriction factor. For a circular pore tube with varying cross-sectional area A(x), the constriction factor is given by thus the constriction factor is related to the variation of the square of the cross-sectional area, and thereby to the square of the body-throat aspect ratio. Despite that the Stokes equations describe the flow during slow displacement, the resulting constriction factor from this slow flow regime is not linearly linked to the aspect ratio suspected to correlate with the slow flow displacement efficiency. A solution to the Laplace equation, ∇ 2 f = 0, would yield a flux ∇ f that is proportional to the cross-sectional area for a cylindrical tube. Many physical processes are governed by the Laplace equation, such as electrical conductance and steady-state diffusion. For electrical conductance, the fractional relation between the effective conductivity σ e of a porous medium to the constant conductivity in the pore space σ can be described as where τ c is an electric field line-based tortuosity, while C c is a constriction factor [34]. The inverted fractional relation σ /σ e is called the formation resistivity factor [35]. For a circular pore tube the constriction factor is thus the electrical constriction factor is related to the variation in cross-sectional area A(x), and thereby to the pore bodythroat aspect ratio. The electrical constriction factor and the displacement efficiency for the primary drainage process is thus related to the same geometrical feature of the pore space, i.e., the aspect ratio between pore throats and pore bodies. We can therefore expect that the electrical constriction factor and the displacement efficiency are related to each other. Consequently we cross-plotted the electrical constriction factor versus the primary drainage displacement efficiency, as shown in Fig. 10. As seen from this plot, there is a clear correlation between the constriction factor and the displacement efficiency. It is therefore postulated that the electrical constriction factor can be used to estimate the displacement efficiency.
For natural porous media, the variation in tortuosity τ c is limited and correlated to porosity, while the porosity φ and formation resistivity factor σ /σ e can be measured by conventional laboratory measurements or estimated from well logs. Thus, we can obtain estimates for the electrical constriction factor C c , both through laboratory measurements and from well logs, and C c can sequentially be used to estimate the displacement efficiency.

VIII. THERMODYNAMICALLY BASED WETTABILITY MEASUREMENTS
In this section we will return to the thermodynamically based estimation for the wettability. During an immiscible displacement the applied work W can be measured externally. Through the calculated displacement efficiency η we then get an estimate for the change in the free energy F through Eq. (12): The change in free energy from the change in surface areas is given as [see Eq. (9)]: Combining Eq. (30) and Eq. (31) we get The fluid-fluid interfacial tension σ can be considered a pure fluid-fluid property which can be measured outside the porous medium. Through pore-scale imaging the surface areas can be considered known, and from external pressure and flow rate measurements also W is known. If we assume that the there is a weak dependency between the efficiency η and the contact angle, then η can be considered a geometrical property obtainable from the images. In this case the only unknown is σ d : This yields a thermodynamically based wettability index as with the contact angle θ obtained from cos(θ ) = ω i . We will now consider an imaged two-phase displacement experiment [33]. In this experiment a sample of sintered soda lime beads are mounted into a core holder, and placed inside a micro-CT for imaging. The experiment started with the sample being completely filled with water, and then slowly displaced by n-dodecane. The external fluid pressures were continuously monitored, and the fluid distribution was imaged at different saturation steps. For a more in-depth description of the experiment, we refer to Ref. [33]. A cross-section of the sintered soda lime beads sample is visualized in Fig. 11. This pore space was also included in the calculations shown in Figs. 9 and 10. In Fig. 11 we have also plotted the inverse of the measured external pressure when the different parts of the pore space was invaded by the injected n-dodecane.
The core holder is shown as a hatched area in the left and right boundary of Fig. 11. One can observe that the displacement pressure is lower along the core holder, and one could argue that the contact angles on the core holder is different from the contact angles on the beads. For our calculations we therefore assumed the beads to be water wet, while the core holder was assumed intermediate wet.
We simulated the displacement similar to the process described in Sec. V. However, we used a contact angle of 90 degrees versus the core holder. The resulting connected parts for different radii values r are plotted in Fig. 12. Note the correlation with the inverse pressures shown in Fig. 11. The efficiency η for the different drainage steps in the experiment was obtained as where s ⊂ is the subspace of the pore space where ndodecane replaces water during saturation step s. The external work is estimated as W s = φ s s w p s , where p s is the externally measured pressure difference at the end of step s, while φ s s w is the change in water volume. The changes in surfaces A 12 and A 2s are also measured from the threedimensional images of the fluid distribution.
We can then calculate the wettability index as given by Eq. (34) for the displacement process. The resulting contact angles θ = cos −1 (w i ) are plotted in Fig. 13. In the same plot we have the resulting contact angles without including the efficiencies η s . Even though this sample has a high efficiency, as seen in Fig. 10, there is still a significant difference between the results when invoking the efficiency and without.
The inclusion of efficiency seems to reduce the spread in contact angles, however, there is a significant spread for the low and high saturation values. For the low values more of the displacement occurs close to the core holder, and might be affected by our assumption of the difference in wettability between the beads and the core holder. For the high saturation values the displacement is short movements of a high number of interfaces, which is expected to be more prone to imaging errors.
The average contact angle is found as 59 degrees, which is just inside the value range 30-60 degrees in reported in Ref. [33], while significantly different from the thermodynamically based wettability calculation in Ref. [9] giving a contact angle of 43 degrees. Excluding the efficiency we obtain a 033113-11 contact angle of 46 degrees. This is close to the thermodynamically derived contact angle of 48 degrees obtained for the same material without taking into account the efficiency [13]. The same study measured the contact angles just before displacement events, which led to averages of 64 and 68 degrees [13], which is close to the contact angles obtained by our method. The difference between the contact angles of 64 and 68 degrees and the thermodynamically derived contact angle of 48 degrees might be assigned to the efficiency, and in this respect match our results quite well.
We now want to estimate the wettability from measures available from traditional experiments. We then use a full primary drainage cycle to approximate the thermodynamically based wettability from Eq. (34) as where η is the drainage efficiency for the total pore space, while A s is the surface area. Note that from imaging of the pore space of a sample we can obtain estimates for both η and A s . We can also use the correspondence shown in Fig. 10 to estimate the displacement efficiency from electrical measurements. The surface area can also be estimated by several different methods, e.g., from adsorption experiments [36][37][38], or from two-or three-dimensional imaging. Last, the externally applied energy W is obtained by integrating the pressure difference curve as Thus, all parameters in the approximation in Eq. (36) can be obtained by traditional experiments. Using this approximation on the drainage experiment in the soda beads we obtain a contact angle of 63 degrees, which is then compared to the value of 59 degrees given by the calculations above. Thus, for our experiment the wettability approximation given by Eq. (36) with input values which are available from traditional experiments gives a fair match with the more complete calculations.

IX. SUMMARY
In this article we have introduced pore-scale definitions for thermodynamically reversible and irreversible displacement processes, called isons and rheons. These definitions are based on the Laplace pressure p over fluid-fluid interfaces with relation to the interface movement, as given by the integral When the integral in Eq. (38) equals the externally applied work, then the displacement process is an ison, while it is a rheon when the integral in Eq. (38) equals the internal rate of energy dissipation. We have investigated the energy efficiency of displacement processes in simplified porous media consisting of sinusoidal tubes, and observed that for a high number of tubes the energy efficiency converges towards a geometrically determined efficiency. This observation together with experimental and simulated results indicate that for pore structures with a high number of pore bodies, the efficiency of primary displacement processes is a function of the pore space geometry.
The geometry of the pore structure is thus linked to a local displacement efficiency factor, which can be integrated to give the primary displacement efficiency for the full porous medium. We have calculated this efficiency for a range of natural porous media. While the fluid flow in the displacement process is described by the Stokes equation, we argue that the displacement efficiency has a cross-property relation with transport processes described by the Laplace equation. This opens up for estimating the primary displacement efficiency from traditional and widely available measurements.
We have used the local displacement efficiency to obtain a thermodynamically based wettability index for an experiment with pore-scale images of the drainage process. The obtained wettability index is in the vicinity of values obtained by other authors. We also showed how our methodology can be used to approximate the wettability index from traditional experiments, namely from the electrical conductance, the capillary pressure curve and surface area approximation. For our experiment, the obtained approximation is similar to the calculated value.

ACKNOWLEDGMENTS
We thank Steffen Schlüter for sharing the sintered soda lime beads data. This research was supported by the Research Council of Norway through its Centers of Excellence funding scheme, Project No. 262644, PoreLab.

APPENDIX A: DERIVATION OF PORE-SCALE EQUATIONS FOR ISONS AND RHEONS
In this Appendix we will derive pore-scale equations used for describing the thermodynamically reversible and 033113-12 irreversible displacement processes, i.e., the ison and rheon processes. At the interstitial scale, i.e., inside the pore space, a slow (creeping) flow u is assumed governed by the Stokes equation supplemented by the continuity equation: where u is the interstitial (pore-scale) fluid velocity, p is the pressure, and μ is the viscosity of the fluid [22,39]. In the following we will refer to Eqs. (A1) simply as the Stokes equations. We assume a no-slip condition at the fluid-solid interface.
As the pressure is discontinuous at the fluid-fluid interfaces, the Stokes equations are not well defined at these interfaces. Further, movement of fluid-fluid interfaces lead to movement of three-phase contact lines. This will give nonzero fluid velocity at the solid surface, contradicting our no-slip condition. This conflict between the no-slip condition and the moving three-phase contact line leads to singularities in the solution of the Stokes equation [40]. This classical problem is still not fully solved [41]. We will not delve into this problem; we assume that the no-slip condition is relaxed in the vicinity of the three-phase contact line. For all internal points in the two fluids the pressure and thereby the Stokes equations are well defined. As we want the continuity equation to hold universally, we can define the velocity of both fluids equal at the interface. Under this assumption the velocity field u will be smooth, also over the interfaces. Due to possible differences in viscosity, the gradient of the pressure field is not necessarily smooth.
Applying the divergence theorem, and invoking that the fluid velocity u is a solenoidal vector field from the continuity equation Eq. (A1b), we obtain where n is the outward pointing unit normal field of the boundary δ . The boundary δ = δ e ∪ δ s ∪ δ f can be decomposed into the external boundaries of the pore space Here δ i indicates that we are integrating over the boundary for fluid i, thus on only one side of the fluid-fluid interface. Surface properties can be considered scale dependent, e.g., surface roughness can be considered as an effective property representing small scale surface geometry, similarly small scale wettability variation can be treated as an effective surface property. We are only considering heterogeneities at a scale larger than the range of the van der Waals forces, still smaller than what is resolved by our pore structure description ⊂ V . Surface roughness and wettability heterogeneity can then be represented by contact angle hysteresis [42,43]. With such effective representation of smaller scale heterogeneities, the increase in free energy for the fluid-fluid surfaces will have a nonreversible component as given by the hysteresis model. Both surface roughness and varying wettability induce small and fast release of energy from the surfaces, thus violating the assumption of increase in free energy for all interfaces used for the pore-scale definition of an ison process [Eq. (3)]. Such release of energy will increase local flow rates, and thereby contradict that pQ ∇ p · udV . Additionally, incompressibility would be a 033113-13 FIG. 15. The connected part C[X (r)] of the subsets X (r) shown in Fig. 14, equivalently c r p (x) as given by Eq. (23). By connected part, we mean the part that is connected to the top of this twodimensional model. too simple fluid description for representing these rapid irreversible processes. In our derivations we considered these effects of second order, i.e., the process is considered close to reversible.
One could construct a pore space geometry where the free energy released from one fluid-fluid interface was exactly balanced by the increase in free energy of all the other interfaces. This would lead to a thermodynamic reversible internal redistribution of fluids, contradicting that this is a rheon displacement. For such a reversible process we can no longer assume a limited time horizon for the internal redistribution of fluids, thus the external work on the system during the internal displacement is no longer negligible, and thereby violating the assumption pQ ∇ p · udV used in Eq. (4). This is only of theoretical interest, in practice the release of energy from moving a fluid-fluid interface would lead to a thermodynamically irreversible process, and the release of energy from this process would be fast compared to the increase in energy from the externally applied work.