Oil Fragmentation, Interfacial Surface Transport and Flow Structure Maps for TwoPhase Flow in Model Pore Networks. Predictions Based on Extensive, DeProF Model Simulations
Department of Civil, Survey and GIS Engineering, TEI Athens,
Ag. Spyridonos,
Athens
12210  Greece
^{*} Corresponding author email: marval@teiath.gr
Received:
16
November
2016
Accepted:
25
September
2017
In general, macroscopic twophase flows in porous media form mixtures of connected and disconnectedoil flows. The latter are classified as oil ganglion dynamics and drop traffic flow, depending on the characteristic size of the constituent fluidic elements of the nonwetting phase, namely, ganglia and droplets. These flow modes have been systematically observed during flow within model pore networks as well as real porous media. Depending on the flow conditions and on the physicochemical, size and network configuration of the system (fluids and porous medium), these flow modes occupy different volume fractions of the pore network. Extensive simulations implementing the DeProF mechanistic model for steadystate, onedimensional, immiscible twophase flow in typical 3D model pore networks have been carried out to derive maps describing the dependence of the flow structure on capillary number, Ca, and flow rate ratio, r. The model is based on the concept of decomposition into prototype flows. Implementation of the DeProF algorithm, predicts key bulk and interfacial physical quantities, fully describing the interstitial flow structure: ganglion size and ganglion velocity distributions, fractions of mobilized/stranded oil, specific surface area of oil/water interfaces, velocity and volume fractions of mobilized and stranded interfaces, oil fragmentation, etc. The simulations span 5 orders of magnitude in Ca and r. Systems with various viscosity ratios and intermediate wettability have been examined. Flow of the nonwetting phase in disconnected form is significant and in certain cases of flow conditions the dominant flow mode. Systematic flow structure mutations with changing flow conditions have been identified. Some of them surfaceup on the macroscopic scale and can be measured e.g. the reduced pressure gradient. Other remain in latency within the interstitial flow structure e.g. the volume fractions of − or fractional flows of oil through − connecteddisconnected flows. Deeper within the disconnectedoil flow, the mutations between ganglion dynamics and drop traffic flow prevail. Mutations shift and/or become pronounced with viscosity disparity. They are more evident over variables describing the interstitial transport properties of process than variables describing volume fractions. Τhis characteristic behavior is attributed to the interstitial balance between capillarity and bulk viscosity.
© M.S. Valavanides, published by IFP Energies nouvelles, 2018
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Notation
∼: A tilde embellishment indicates a dimensional variable
No tilde embellishment denotes a dimensionless variable
Φ^{′}: A prime indicates a physically admissible value of variable Φ (see Sect. 1.5)
In the following lists the far right columns indicate:
Scale: M (macroscopic), S (statistical average), P (pore scale) (see Sect. 1.2, end)
Definition: corresponding equation (number) or paragraph number
Latin characters
_{Notation}  Interpretation of variable  Scale  Definition 

Porous medium crosssectional area of REV  M  (2)  
Total o/w surface area  S  (43)  
Total o/w surface area in the DTF domain  S  (41)  
Total o/w surface area in the DOF domain  S  (42)  
Total o/w surface area of ganglia  S  (36)  
Oil/water surface area of mobilized drops, DTF domain  S  (39)  
Oil/water surface area of mobilized ganglia  S  (33)  
Mobilized o/w surface area of mobilized ganglia  S  (31)  
Oil/water surface area of stranded drops, DTF domain  S  (40)  
Oil/water surface area of stranded ganglia  S  (34)  
Stranded o/w surface area of mobilized ganglia  S  (32)  
,  Mean (across all unit cells) surface area of oil/water interfaces, advancing and receding  P  (27) 
_{}  Mean (across all unit cells) surface area of o/w interfaces (mobilized and stranded) of any iclass mobilized ganglion  S  (29) 
Mean (across all unit cells) surface area of the mobile oil/water interface of any mobilized ganglion  S  (27)  
_{}  Mean (across all unit cells) surface area of the oil/water interface of any iclass stranded ganglion  S  (30) 
_{}  Mean (across all unit cells) surface area of the stranded oil/water interface of any iclass mobilized ganglion  S  (28) 
Ca  Capillary number  M  (1) 
Fraction of classi core (‘C’) ganglion cells − conducting  S  (20)  
Fraction of classi nonconducting (‘X’) ganglion cells  S  (21)  
Fraction of classi end (‘E’) ganglion unit cells − conducting  S  (22)  
D_{Cj}  Diameter of jclass chamber  P  2.1 
D_{Ti}  Diameter of iclass throat  P  2.1 
f_{EU}  Energy efficiency of the process  M  (58) 
f_{j}  Occurrence probability of any iclass network element  P  2.1 
_{}  Coefficient of oil fragmentation  S  (48) 
,  Conductance of the various unit cells in the GD and DTF domains  P  1.3 
I_{max}  Maximum ganglion size class  P  1.4 
Absolute permeability of the pore network  M  (5)  
ℓ  Network lattice constant  P  2.1 
m^{D}  Mean mobilization probability of drops  S  2.3 
Mean mobilization probability of iclass ganglia  S  2.3  
Number density of iclass ganglia  S  (15)  
N_{d}  Network dimensionality  P  2.1 
Number of pore unit cells occupied by each iclass ganglion  P  (18)  
Average number of oil drops p.u.v. of DTF domain  S  (38)  
Critical number density of minimum size oil droplets  P  (46)  
,  Occurrence probabilities of unit cells in the GD and DTF domains  S  (13) 
,  Macroscopic flow rates of oil and water  M  (2) 
Mean flow rate of oil through CPF  M  (49)  
Flow rate within any conducting iclass ganglion unit cell  P  (19)  
r  Oil/water flow rate ratio  M  (2) 
S_{w}  Mean water saturation (entire network)  M  1.2 
S_{o,D}  Mean oil saturation in the DTF domain  S  (9) 
S_{o,G}  Mean oil saturation in the GD domain  S  (9) 
Mean ganglion velocity of iclass ganglia  S  (12)  
Superficial velocity of oil in the DTF domain  M  (50)  
Superficial velocity of oil in the DOF domain  S  (12)  
Superficial velocity of oil in the GD domain  S  (50)  
Superficial velocity of oil in the CPF  M  (6)  
,  Superficial velocity of oil and water  M  (2) 
Superficial velocity of the o/w interface in the GD domain  S  (52)  
Superficial velocity of the o/w surface in the DTF domain  S  (53)  
Superficial velocity of the o/w interface in the DOF domain  M  (54)  
U_{w}  Superficial velocity of water  M  (2) 
U_{w,DOF}  Superficial velocity of water in the DOF domain  M  (7) 
Pore network reference volume  M  1.2  
V_{D}  Average (mean) oil drop volume  S  (37) 
_{VD,1}  Minimum size of an oildroplet  P  (44) 
Volume of iclass ganglia pore unit cells  P  (24)  
Mean size of all ganglia (mobilized and stranded)  S  (24)  
Mean size of the mobilized ganglia  S  (25)  
Mean size of the stranded ganglia  S  (26)  
Size of the ensemble of physically admissible solutions  M  (16)  
W  Rate of mechanical energy dissipation within the process  M  (56) 
_{}  Mechanical power dissipation for an equivalent 1phase flow  M  (56) 
W_{C}  Mechanical power dissipated in CPF  M  (57) 
W_{D}  Mechanical power dissipated in DTF  M  (57) 
W_{G}  Mechanical power dissipated in GD  M  (57) 
x  Macroscopic pressure gradient  M  (5) 
x_{pm}  Set composed of pore network parameters  _{P}  _{1.1} 
Direction of macroscopic flow  M  (5) 
Greek characters
Notation  Interpretation of variable  Scale  Definition 
β  Porous medium volume fraction occupied by the connected oil  M  (8) 
ζ  Parameter measuring the sharpness of the ganglion size distribution  S  (15) 
Fractional flow through connected pathways  M  (49)  
Fractional flow through drop traffic flow  M  (49)  
Fractional flow through ganglion dynamics  M  (49)  
θ_{A}, θ_{R}  Advancing and receding contact angles of the o/w interface  P  1.1 
κ  Oil/water viscosity ratio  M  (3) 
Oil/water interfacial tension  P  (1)  

Viscosity of water  M  (3) 
Viscosity of oil  M  (3)  
Fractional transport of the o/w interfaces through DTF  S  (55)  
Fractional transport of the o/w interfaces through GD  S  (55)  
σ_{3D}  Mean coordination number of pore network  P  2.1 
χ_{G}  Tortuosity of ganglia in the pore network  P  2.1 
ω  Fraction of all the ganglion cells over all the DOF region cells  M  (9) 
Abbreviations
CPF  Connected Pathway Flow 
DeProF  Decomposition in Prototype Flows 
DOF  Disconnected Oil Flow 
DTF  Drop Traffic Flow 
FAV  Flow Arrangement Variable(s) 
GD  Ganglion Dynamics 
GCD  Ganglion Cells Domain 
o/w  oil/water 
p.m.  porous medium 
p.u.v.  per unit volume 
p.u.v.p.m.  per unit volume of porous medium 
REV  rRepresentative Elementary Volume 
u.c.  unit cell 
Introduction
Fluid flow in porous media, i.e. porous rocks, fracture networks and granular media, is subject to considerable interdisciplinary research activity in physics, earth sciences and engineering. Examples of natural and engineered processes include hydrocarbon recovery, carbon dioxide geosequestration, soil drying/wetting or pollution remediation, operation of packedbed reactors, operation of proton exchange membrane fuel cells, etc.
Immiscible twophase flow in porous media is a process whereby a resident phase is replaced by − or flows simultaneously with − another phase under the effect of externally imposed pressure gradients, induced artificially (e.g. pumping) or naturally (e.g. gravity). One of the fluids is wetting the solid (“water”) in presence of the other, nonwetting, fluid (“oil”). One possible flow mode is the simultaneous, concurrent, forced flow of oil and water at preselected flow rates. Other possible flow modes could be imbibition, that is, displacement of oil by water, and drainage, that is, displacement of water by oil.
These processes take place across widely different length scales. Hydrodynamic flow instabilities and pore scale disorder induce flow disconnection and create a hierarchical system typically resulting in complex flow patterning. The disconnection of the nonwetting phase (oil) produces small fluidic elements, called ganglia and droplets. Ganglia are microfluidic elements (oil blobs of very small size) that occupy the void volume of a few (one to ten) pores. Droplets are even smaller fluidic elements, so that any pore may contain many of these. When combined, the physicochemical properties of the wetting, the nonwetting and the solid phase Porous Medium (p.m.) create interfaces separating the nonwetting fluidic elements from the wetting phase. The collective effect of inherent flow mechanisms, prevailing over the pore (microtomillimeter) scale domain and described by Stokes and YoungLaplace equations, results in the peculiar phenomenology of twophase flows in porous media observed on the core, fracture, or field/production scales.
In general, macroscopic twophase flows in porous media form mixtures of connected and disconnectedoil flows. The latter are described as oil Ganglion Dynamics (GD) and Drop Traffic Flow (DTF), modes that have been systematically observed during flow within model pore networks (Avraam and Payatakes, 1995; Tzimas et al., 1997; Tallakstad et al., 2009; Krummel et al., 2013; Datta et al., 2014; Armstrong et al., 2016). Apart from laboratory studies, virtual studies, implementing dynamic pore network simulations (Constantinides and Payatakes, 1996 and 2000; Knudsen et al., 2002; Knudsen and Hansen, 2002; AlGharbi and Blunt, 2005; Nguyen et al., 2006), or lattice Boltzmann methods (Pan et al., 2004; Ghassemi and Pak, 2011; Ramstad et al., 2012; Armstrong et al., 2016), have also addressed disconnectedoil flow modes. With recent advances in imaging technology, these flows are observed in real porous media as well (Georgiadis et al., 2013; Oughanem et al., 2015; Rücker et al., 2015). The observed flow structures show a significant and − presumably − systematic mutation ranging from practically no oilflow, to large and small oil ganglion dynamics, to drop traffic and connected pathway flow mixtures, depending − primary − on the flow conditions and − secondary − on the physicochemical, size and network configuration of the oilwaterp.m. system (Avraam et al., 1994; Avraam and Payatakes, 1995 and 1999; Krummel et al., 2013; Armstrong et al., 2016). These studies show that during simultaneous flow of oil and water, the disconnected phase (oil) contributes significantly to the total flow and, in certain cases even exclusively, e.g. emulsion flow (Cobos et al., 2009; Guillen et al., 2012). Furthermore, the flow rate against pressure gradient relation is found to be strongly nonlinear, and to be strongly affected by the physical parameters that pertain to the fluidfluid interfaces.
Focusing on steadystate, immiscible, twophase flow on porous media, it is regulated by the incessant exchange of position and momentum, between disconnected, nonwetting fluidic elements, as these migrate downstream within the pore network at macroscopically constant flow rates under macroscopically constant pressure gradients. The average flow is a mixture of structurally different flow configurations, or ‘prototype flows’ as identified by Avraam and Payatakes (1995), namely ConnectedOil Pathway Flow (CPF) and DisconnectedOil Flow (DOF), comprising GD and DTF. In the same work, a map of the flow regimes − according to the prevalent prototype flow − has been drafted to coarsely demarcate the dependence of the flow structure on the flow conditions.
To elucidate the actual, inherent flow mechanisms that form the structure of flow regimes and regulate the overall (macroscopic) behavior of immiscible steadystate twophase flows in pore networks, Valavanides and Payatakes (2000 and 2001; also Valavanides, 1998 and 2012) developed the DeProF mechanistic model. The model is based on the concept of decomposition in prototype flows (see next section). It accounts the porescale mechanisms and the networkwide cooperative effects. The sources of nonlinearity, mainly caused by the motion of populations of oil/water interfaces, are modeled satisfactorily and their effect is projected in the macroscopic scale by combining effective medium theory with appropriate expressions for poretomacro scale consistency of the mass transport of individual phases. DeProF detects all the physically admissible, interstitial flow configurations. The latter form an ergodic ensemble defining the average configuration of the macroscopic flow. Reduced pressure gradient and relative permeabilities are estimated as ensemble averages. DeProF can be efficiently used as a simulation tool. The quantitative and qualitative agreement between existing sets of data and the corresponding theoretical predictions of the DeProF model is excellent (Fig. 3 in Valavanides and Payatakes, 2000; Fig. 5 in Valavanides, 2012). As a simulation tool, it has revealed latent critical characteristics of the flow structure. Its predictions on energy efficiency aspects have been recently validated by a major retrospective examination of independent laboratory studies (Valavanides et al., 2016).
Using DeProF, it is possible to predict, appropriately defined key physical quantities to describe the interstitial structure of the flow on different levels of specificity: ganglion size and ganglion velocity distribution, mean ganglion size, fractions of mobilized/stranded populations, specific (per unit volume of p.m.) surface, velocity and volume fractions of mobilized and stranded interfaces, oil fragmentation, etc.
The scope of the present work is to elucidate the inherent phenomenology of twophase flow in pore networks, to understand the role of disconnected flow and to reveal if any dynamic similarities apply across different oilwaterp.m. systems. To this end, extensive simulations implementing the DeProF algorithm for systems similar to those studied by Avraam and Payatakes (1995; 1999) have been carriedout to derive maps that describe the dependence of the flow structure on the independent flow variables, namely the capillary number, Ca, and the flow rate ratio, r, or, equivalently, the reduced, superficial velocities of oil and water. The simulations span 5 orders of magnitude in Ca(− 8 ≤ logCa < − 3) and r(− 2 ≤ logr ≤ 2). Systems with various − favorable and unfavorable − viscosity ratios have been examined. The novelty of the present work is the delivery of maps of physical quantities that describe the interstitial flow structure of twophase flows in a typical 3D model pore network of the chamberandthroat type, over broad domains of capillary number, flow rate ratio and viscosity ratio values. These are produced by extensive and systematic simulations with the DeProF mechanistic model. Previous works by Valavanides and Payatakes (2000; 2001; 2002; 2003; 2004) addressed the performance of the DeProF model, benchmarked the specificity of its predictions and revealed critical flow characteristics over a single order of magnitude on Ca and for systems with specific viscosity ratio values − to simulate the Avraam and Payatakes (1995) experiments. Using these extensive simulations and the associated flow structure maps, it is now possible to understand the mechanisms underlying the significant mutations observed in twophase flows in porous media in mixtures of connected and disconnected oil flow (connectedoil pathway flow, small and large ganglion dynamics, drop traffic flow) and reveal the systematic nature of these mutations. If such trends can be revealed, then, we would be more confident in associating the macroscopic flow behavior to the interstitial flow structure. The former can be relatively easily observed/measured while the latter requires sophisticated equipment and/or elaborate methodologies to be evaluated.
Why do we need to understand the behavior of the flow in terms of its internal structure? Because disconnection of the nonwetting phase (oil) is a critical phenomenon; it enriches the flow with many internal degrees of freedom (by expanding the ensemble of physically admissible flow configurations) and creates an abundance of interface surfaces between wetting and nonwetting phases. The latter have a strong impact on the momentum balance as they act as buffers and/or as instigators of local hydraulic perturbations, with mobilization of the nonwetting phase − an important characteristic in enhanced oil recovery processes (Alvarado and Manrique, 2010). In addition, interface surfaces play an important role as many physicochemical processes with industrial relevance take place on or across interfaces, e.g. the preferential adsorption at the o/w interface of tracer nanoparticles (Ryoo et al., 2010) and other species (SerresPiole et al., 2012), or the rate of dissolution of the nonwetting in the wetting phase (Sahloul et al., 2002) directly related to the surface area per unit p.m. volume, have a strong impact on the spreading of contaminants in groundwater aquifers, the storage of CO_{2} in brine filled formations, the efficient operation of PEM fuel cells (Bazylak et al., 2008) or tricklebed reactors (Van de Merwe and Nicol, 2009) and other industrial applications.
The paper deploys in four main sections. Section 1 provides an overview of the essentials of the DeProF model. Section 2 describes the basic characteristics of the model pore network and the associated modeling provisions with regards to the discrete fluidic elements (ganglia and droplets) that comprise the DOF. Section 3 lists the key variables that have been defined to describe the interstitial structure of the flow, e.g. specific surface area and surface transport of the o/w interfaces, oil flow partitioning (fractional flows) through prototype flows, ganglion size distribution, extent of oil disconnection, etc. The reader who is familiar with the DeProF model may skip Section 2 (or even Sect. 3) and proceed directly to Section 4 that presents and discuss the results of the simulations. Finally, general conclusions are furnished as a last section.
1 The DeProF Hierarchical mechanistic model for steadystate twophase flow in pore networks
The mechanistic model DeProF developed by Valavanides and Payatakes (2000, 2001, 2012), can be used to obtain the solution to the problem of steadystate twophase flow in pore networks, by predicting the reduced macroscopic pressure gradient, x, given the values of physical quantities defining the flow conditions and the oilwaterp.m. system properties.
1.1 The independent variables for steadystate twophase flow in pore networks
The conditions of steadystate twophase flow in pore networks are defined by the values of the flow rate intensities, i.e. superficial velocities, and , conventionally defined as oil and water flow rates, and , divided by the p.m. crosssectional area, . The set of superficial velocities may be appropriately reduced and replaced by a set of dimensionless variables, comprising the capillary number, (1) and the oil/water flow rate ratio, (2)where is the viscosity of water and is the oil/water interfacial tension.
The oilwaterp.m. system is characterized by a set of physical parameters describing key physicochemical and structural (network geometry, topology) properties. This set comprises:

the absolute permeability of the p.m., _{;}

the oil/water viscosity ratio;

the dynamic advancing and receding contact angles of the o/w interface on the pore walls, θ_{A} and θ_{R}, and;

a (sub)set, x_{pm}, comprising all the dimensionless geometrical and topological parameters of the pore network that affect the flow (e.g. porosity, genus, coordination number, normalized chamber and throat size distributions, chambertothroat size correlation factors, etc).
1.2 The concept of Decomposition into Prototype Flows (DeProF)
In the general case of steadystate, onedimensional, twophase flow in pore networks, depicted by the left sketch of Figure 1a, oil and water are continuously supplied in the p.m. along the macroscopic flow direction, , with constant flow rates, and , so that the dimensionless parameters defining the flow conditions, Ca and r, attain constant values. This could be the case of using syringe pumps for the continuous, simultaneous displacement of oil and water through a core.
Figure 1 (a) Schematic representation of the actual flow (left sketch) and its theoretical decomposition into prototype flows: Connectedoil Pathway Flow (CPF) and DisconnectedOil Flow (DOF) (right sketch), (b) A microscopic scale representation (snapshot) of a DOF region. An oil ganglion of size class 5 is shown. For simplicity, all cells are shown identical and the lattice constant is shown expanded (chambers and throats have prescribed size distributions). The dashed, light grid lines define the network unitcells. The thick dashed line separates the Ganglion Dynamics (GD) cells domain and the DropTraffic Flow (DTF) cells domain. 
Water is considered to be the wetting phase and always retains its connectivity, while oil (as the nonwetting phase) gets disconnected to form fluidic elements of various sizes: from infinitelong and slim connectedoil pathways, to short ganglia, to tiny oil droplets. The degree of wettability may vary and this can be accounted for by the magnitude of the contact angles of the oil/water interfaces on the pore walls. The average water saturation is S_{w}.
The macroscopic flow can be decomposed into two − structurally different − flow patterns,

Connectedoil Pathway Flow (CPF) and;

Disconnected Oil Flow (DOF).
The latter can be further decomposed in

Ganglion Dynamics (GD) flow and;

Drop Traffic Flow (DTF).
All types of flow have been observed experimentally in model porous media (Avraam and Payatakes, 1995, 1999; Tsakiroglou et al., 2007; Tallakstad et al., 2009; Krummel et al., 2013; Armstrong et al., 2016) as well as in real porous media (Van de Merwe and Nicol, 2009; Georgiadis et al., 2013; Oughanem et al., 2015; Rücker et al., 2015).
In the DeProF model, each prototype flow (CPF, GD and DTF) has the essential characteristics of the corresponding flow pattern in suitably idealized form and so, the pore scale mechanisms are incorporated in the prototype flows. In the Connected Pathway Flow (CPF) region the oil retains its connectivity and flows with virtually onephase flow. The p.m. volume fraction occupied by the connected oil is denoted by β. The DisconnectedOil Flow (DOF) regime is defined as the region composed of the rest of the pore network, so the associated volume fraction equals (1 − β). This modeling decomposition is represented schematically by the right sketch in Figure 1a, where the volumetric fractions of the prototype flows have been isolated and reaccounted over separate parts of the pore network Representative Elementary Volume (REV) of size . The size (volume), , of the REV is selected so that, starting from a critical (minimum /threshold) size, any measurable quantity, e.g. pressure differences and flow rates of oil and water, interstitial structure of flow, do not change if the REV is further increased. The problem with selecting the size of the REV is identical with the problem of selecting a representative core size of a real p.m. Note however that, due to the inherent degrees of freedom instigated by the disconnection of the nonwetting phase, the REV pertaining to immiscible twophase flow processes should be much larger that the REV pertaining to saturated flow − even if the same pore network is considered.
A microscopic scale representation (a snapshot) of a typical DOF region is shown in Figure 1b, where an oil ganglion having a typical ‘cruising’ configuration (Valavanides and Payatakes, 2001) is shown at the center. All the cells that accommodate parts of this (or any other) oil ganglion are called ganglion cells and are demarcated with a thick dashed line. The rest of the cells in the DOF regions are cells containing water and oil drops. These cells comprise the regions of the GD and DTF domains respectively.
The fraction of all the ganglion cells over all the DOF region cells is denoted by ω, and is called the GD domain fraction. The DTF domain fraction in the DOF region equals (1 − ω).
S_{w}, β and ω are called Flow Arrangement Variables (FAV) because they provide a coarse indication of the structure (or ‘immiscible mixing’ of the constituent prototype flows) of the flow pattern. One of the objectives of the DeProF model algorithm is to determine the values of S_{w}, β and ω that conform with the externally imposed flow conditions (Ca, r).
The flow analysis is carriedout on two length scales, the macroscopic scale (10^{12} pores, or more, the REV size) and the microscopic (or pore) scale, extending within a single or over a few pore unit cells. DeProF is essentially a stochastic model and the flow analysis produces a system of equations that includes macroscopic water and oil mass and momentum balances and microtomacroscopic scale consistency relationssee next subsection.
Scalewise, the variables used in the DeProF model can be classified into three groups.

Group (P) variables defined at the scale of one or a few pore unit cell(s). These may describe either pore network geometrical characteristics (e.g. any pore chamber diameter), or physical quantities that depend entirely on the pore characteristics (e.g. the surface area or curvature of an o/w interface − receding or advancing − within a unit cell);

Group (S) variables pertaining to Statistical average (scaledup) values of porescale (P) variables, (e.g. the statistical average velocity of an iclass ganglion as it migrates downstream all possible combinations of pore unit cell sizes);

Group (M) variables pertaining to the Macroscopic scale (M) or REV scale. Group (M) comprise two types of variables depending on their principal definition: (a) variables defined/measured directly at the macroscopic scale, e.g. the water saturation, the pressure gradient, the volume fraction of connected oil pathways, the superficial velocities of oil and water, etc. (b) Group (S) variables.
All variables used in the manuscript are extensively listed in the NOTATION.
1.3 The DeProF model equations
Using the DeProF model, one can obtain the solution to the problem of steadystate twophase flow in porous media in the form of the following transfer function (4) where, (5)is the reduced macroscopic pressure gradient, i.e. the actual pressure gradient divided by the pressure gradient of an equivalent onephase flow of water at superficial velocity equal to [the 2nd component of the product in Equation (5)].
Referring to Figure 1, the externally imposed flow rates of oil, , and water, , are partitioned in the prototype flows. In particular, , are the mean oil flow rates in CPF and DOF and the mean water flow rate in DOF (as CPF is oilsaturated, ). The following physical quantities are also assigned to the prototype flows: the mean oil saturations S_{o,G}, S_{o,D}, in GD and DTF respectively, the mean water saturation, S_{w,DOF}, in DOF, the reduced oil mean superficial velocity in CPF, defined by , and the reduced oil and water mean superficial velocities defined by , where m = o, w.
The flow analysis produces the following system of equations (Valavanides and Payatakes, 2000; 2001) (6) (7) (8) (9) (10) (11) (12) (13) where _{} is the reduced number density of the iclass ganglia and I_{max} is the upper truncation ganglion size class. _{} is the number of unit cells occupied by any iclass ganglion. _{} is the reduced mean ganglion velocity iclass ganglia. _{} is the mean mobilization probability of the iclass ganglia. _{} and are the occurrence probabilities of the various unit cells in the GD and DTF domains respectively. _{} and are the reduced conductance of the various unit cells in the GD and DTF domains respectively.
Note that m^{D}, , , and depend on x, as described by semianalytical expressions of momentum balance derived over the pore scale. , , and also depend on the reduced pressure gradient, x, as well as on the unit cell geometry (index ‘a’), the twophase flow pattern in it (index ‘b’) and the ganglion size (index ‘i’). and depend on x and the unit cell geometry (Valavanides, 1998; Valavanides and Payatakes, 2001).
Equation (6) describes the flow within the connected oil pathways (CPF). Pressure gradient and flow rate intensity are essentially related with a Darcytype expression within the oil saturated CPF fraction of the pore network. The next 3 equations describe oil and water balances accounted at the pore scale and at the macroscopic scale, within the REV, ; in particular, the water flow rate balance, Equation (7), the oil flow rate balance, Equation (8), and the oil and water volumetric balance, Equation (9). Equations (10)–(12) express the consistency in accounting the volumes and flow rates of oil and water over the pore and macroscopic scales; in particular, the consistency between the number of gangliaoccupied pore unit cells and the GD volume fraction, Equation (10), the oil volume (mass) as contained in the ganglia distribution and as saturating the GD domain, Equation (11), and the collective transport of migrating swarm of oil ganglia and the flow rate in DOF, Equation (12). In fact these last three equations instigate an abundance of degrees of freedom by allowing the process to attain various ganglion size distributions and deliver a variety of physically admissible solutions − see Section 1.5. Equation (13) results from application of Effective Medium Theory (Kirkpatrick, 1973) to the equivalent onephase flow in the DOF (GD and DTF) region and implicitly represents the transfer function for this region.
The above system of 8 equations relates (I_{max} + 6) dependent variables, namely x, U_{o,CPF}, U_{o,DOF}, U_{W,DOF}, S_{o,G}, S_{o,D} and , along with the three FAV, namely S_{w}, β, ω. The system is closed by imposing an appropriate type of (sharply decreasing) distribution function for the ganglion volumes, Equation (15), that is based on modeling assumptions presented in the following paragraph. This way the system reduces to a single algebraic equation, as Equations (6)–(12) are all semianalytically and/or numerically infused in (13).
Note that the most difficult part in implementing the DeProF model is in estimating the dependence of , , and on the given pore network geometry. This may be accomplished either by furnishing appropriate analytical expressions (if possible) for , , and or by delivering appropriate maps implementing advanced numerical techniques (CFD or LB) and/or sophisticated laboratory experiments. Analytical expressions have already been derived in the past for the structure/geometry of the chamberandthroat type model pore network described in Section 2 (Valavanides, 1998; Valavanides and Payatakes, 2001). These have been used in the present work.
1.4 The Discrete Fluidic Elements of the Non Wetting Phase
Discrete fluidic elements are formed by the disconnection and immiscible dispersion of the nonwetting phase (oil) within the wetting phase (water), as a result of the combined interaction of hydrodynamic and capillary forces during the tortuous flow path within the pore network conduits. The extent of disconnection (or fragmentation) of the nonwetting phase depends on the imposed flow conditions and the physicochemical properties of the two phases (fluids) in the pore network. The disconnected fluidic elements comprise oil blobs of different sizes − separated by continuous water. Their size extends from one to infinitemany pore unitcells.
In general, three classes of discrete fluidic elements are considered, depending on their size: connectedoil pathways comprising the Connected Pathway prototype Flow (CPF), oil ganglia, occupying onetomany unitcells comprising the GD prototype flow, and oil droplets dispersed within water saturated unitcells comprising the Drop Traffic prototype Flow (DTF). As explained in the previous Section 1.2, the three prototype flows occupy the p.m. space in corresponding volume fractions.
The CPF comprises thin, indistinguishable ConnectedOil Pathways (COP), i.e. connectedoil unitcells with infinite length. COPs may occupy or arrange themselves within any β region of the pore network and remain aligned to the macroscopic flow. Of a total number of M unitcells in a reference volume, βM unitcells are occupied by connected oil (the ensemble of ConnectedOil Pathways) and allow the rest (1 − β) M of the unitcells to host the DisconnectedOil Flow (DOF) unitcells, Figure 1a.
DOF occupies the remaining unitcells within any (1 − β)region of the pore network. The DOF domain unit cells are partitioned into two subpopulations: Unit cells that are occupied by a mixture of ganglia of particular size distribution and unit cells containing water oildroplets and water. Oil ganglia unitcells take up a volume fraction of (1 − β) ω of the pore network, leaving to the DTF cells a complementary volume fraction of (1 − β) (1 − ω). The aforementioned partitioning is described by Equations (7)–(9). The DTF unit cells are all considered to be indistinguishable − refer to Section 2.2 and Figure 7.
Ganglion unitcells [Figs. 1b and 2] are not indistinguishable in the sense that each ganglion unitcell is part of a certain ganglion size class. Oil ganglia have a sharply decreasing volume size distribution, defined as (14) where N_{d} is the network dimensionality (=3 for 3D networks), is the average volume of unit cells and denotes the ratio of the total number of iclass ganglion cells over the total number of all ganglion cells in the DOF region. The distribution of the population density of ganglia of size i, is considered to be given by (15)where n_{1}, n_{2}, I_{max} and ζ are parametric values that define the particular physically admissible flow configuration. I_{max} is the maximum attainable size of a ganglion and 0 < ζ < 1 provides a measure of the sharpness of the ganglion size distribution, see Figure 2.
Figure 2 The model of the ganglion size distribution expressed through equation (15) by three parameters, 
The type of ganglion size distribution is dictated by the physics of GD, and is in compliance with experimental observations (Avraam and Payatakes, 1995; Tallakstad et al., 2009) and numerical /pore network simulations (Valavanides et al., 1998). In reality, large ganglia cannot survive because, as they migrate downstream the tortuous paths of the pore network, they breakup into smaller ones. The values of the parameter used in Equation (15) can be estimated either numerically, implementing the DeProF model algorithm or in the laboratory, implementing any of the modern imaging/scanning technologies (Fusseis et al., 20142; Youssef et al., 2014; Georgiadis et al., 2013; Datta et al., 2014; Oughanem et al., 2015; Armstrong et al., 2016).
Indicative DeProF model predictions of physically admissible ganglion size distributions, , are presented in the diagrams in Section 3.1, along with the corresponding values (tabulated) of the FAV (S_{w}, β, ω). The expected, average values of the flow arrangement variables, (S_{w}, β, ω), for different values of the capillary number, the flow rate ratio and the viscosity ratio, are presented in the diagrams in Section 3.2.
1.5 Physically Admissible Flow Configurations
A core feature of the DeProF model is the detection of all flow configurations that are physically admissible under the imposed macroscopic flow conditions. The actual flow at any REV of the p.m. (extending at an order of 10^{9} pores along the main flow direction), ‘wanders’ within the domain of physically admissible flow configurations, ‘visiting’ any one with equal probability or frequency. This equiprobability is a fundamental assumption of the DeProF model, related to the ergodic nature of the process (see ‘2.3. Discussion on Assumption of the Process Ergodicity’ in Valavanides and Daras, (2016)
To detect these flow configurations, the domain {S_{w}, β, ω} of all possible flow arrangement values is partitioned using sufficiently fine steps to obtain a 3D grid (Fig. 3). Then, for each combination of the flow arrangement values (S_{w}, β, ω), the DeProF algorithm is asked to detect (if) a solution of the DeProF Equations (6)–(13) and (15) (exists). In such case the selected triple of {S_{w}, β, ω} values, is allowed as a physically admissible solution (PAS) and is denoted as . Otherwise, the triple is rejected. At the end of the discretization and scanning process, the subdomain of the {S_{w}, β, ω} space, denoted Ω_{PAS}, is detected as the set of physically admissible solutions. This set is depicted by the cloud of small red balls in the diagrams of Figure 3.
Figure 3 Typical domains of physically admissible flow configurations, as predicted by DeProF model simulations. For any fixed, externally imposed, flow conditions, (Ca, r), the twophase flow visits a canonical ensemble of physically admissible flow configurations depicted by the cloud of small, red balls. Here, Ca = 1.19 × 10^{−6}, κ = 1.45 and flow rate ratio values, r, span 2 orders of magnitude. A unique triplet of values (S_{w}, β, ω), the large, black ball, is obtained by averaging over the ensemble and defines the mass center of the red cloud. Note that the shape and position of the red cloud change with flow conditions. 
The measure (“volume”) of Ω_{PAS} is given by (16) where Ω_{PAS}(Ca, r) stands for integration carried over the physically admissible ranges in {S_{w}, β, ω} for the imposed pair of (Ca, r) values. The volume of Ω_{PAS} is a measure of the degrees of freedom of the process and the variety of microstates it can attain; it is also related to the rate of entropy production (configurational entropy) within the process (Valavanides and Daras, 2016). The volume of the ensemble of interstitial flows is a measure of the internal degrees of freedom of the process; its size and shape depend on the externally imposed macroscopic flow conditions. There are flow conditions for which there are many physically admissible interstitial flow arrangements and there are flow conditions for which the macroscopic flow is confined within a tight ensemble: note the change in position, size and shape of the red cloud with changing flow conditions, Figure 3.
By assuming that each physically admissible solution is visited with the same probability, or frequency (assumption of ergodicity), and averaging over their domain, Ω_{PAS}, a unique solution for the macroscopic flow is obtained. For any quantity, Φ ', the corresponding expected mean macroscopic flow quantity, Φ is defined as (17)
A prime is used to denote physically admissible values of any quantity. The symbol without a prime is reserved for the expected, ensembleaveraged value of that quantity.
As with any macroscopic physical quantity, unique triplet of values for (S_{w}, β, ω) can be obtained by averaging over the domain of physically admissible solutions, using Equation (17), with Φ replaced by S_{w}, β and ω. These values define the ensemble average flow configuration of the macroscopic flow. They are indicated by the large black balls in the diagrams of Figure 4.
Figure 4 Input/output block diagram in the DeProF model implementation. The variables describing the structure of the interstitial flow are presented in Section 3. 
In computational practice, the size of the domain of physically admissible solutions, , as well as any expected mean macroscopic flow quantity, Φ, can be determined by any coarse graining (discretization) procedure used to detect the ensemble of physically admissible solutions and delineate its envelope. The phase space {S_{w} × β × ω} pertaining to the FAV, is partitioned in corresponding N_{S} × N_{β} × N_{ω} phase space voxels. The finer the graining, the smaller becomes the size of these voxels. Nevertheless, the smallest discretization size of voxels has a limit; it cannot become smaller than a minimum, critical size in each prototype flow. Suppose that we detect two solutions and associated with two adjacent voxels, i and i+1. These solutions are virtual in the sense that they have been detected by a numerical algorithm implementing a particular discretization. In essence these have been pickedup from the complete ensemble of the physically admissible solutions. So, there is a great chance that other flow configurations exist − other than those detected (as solutions) by the discretization scheme. This possibility could be verified by rediscretizing the phase space into a finer grid implementing smaller voxels. The new virtual ensemble of physically admissible solutions would eventually contain more solutions. Nevertheless, it would not occupy a volume in the FAV phase space larger than the previous ensemble (associated to the coarse grid). The only difference between the two virtual ensembles would be on the texture of their envelopes. The latter would have a finer surface than the former. We would be able to repeat the rediscretization with finer and finer grids, resulting in smaller and smaller voxels. But we would eventually reach a critical size whereby solutions would not have discernible characteristics. This limiting critical size (threshold) depends on the structure of the pore network (its fractal characteristics) and on the modeling of the prototype flows (the characteristic sizes of the disconnected fluidic elements).
To conclude the present chapter, the DeProF model implementation in terms of input/output is outlined in Figure 5. Imposed flow conditions are provided by three groups of dimensional variables that describe the pore network structure, the physicochemical characteristics of the two immiscible fluids and the solid phase of the p.m. as well as the flow rates with which oil and water are driven through the core. These groups are transformed into two sets of dimensionless variables describing the system parameters and the flow conditions. The latter provide the input to the DeProF algorithm (Sections 1.3–1.5). The output comprises the rheological description of the process, i.e. the reduced pressure gradient in terms of the flow conditions and the system parameters, as well as the structure of the interstitial flow that can be determined by appropriate variables (see Section 3).
Figure 5 Typical jikclass unit cell of the model pore network used in the simulations. The thick inclined arrow represents the orientation of the lattice skeleton relative to the macroscopic pressure gradient. 
Here, a remark must be made with reference to the conventional approach /practice in using the saturation, S_{w}, as the independent variable of the process. It is based on the perception that disconnected oil blobs and other fluidic elements (ganglia and/or droplets) do not move with the average flow but remain stranded in the pore medium matrix. This situation arises when flow conditions of relatively “small values” of the capillary number are maintained. In that case, and that case only, does a linear relation applies between pressure gradient and superficial velocities − an extended Darcy law − whereby the constant of linearity, i.e. the relative permeability, depends on the saturation. Nevertheless, there is evidence that disconnected flow is a substantial and sometimes significant − if not prevailing − flow pattern (see introductory notes), the causeeffect relation describing the flow is strongly not linear and relative permeabilities depend on the flow rates (Knudsen and Hansen, 2002; Nguyen et al., 2006; Tallakstad et al., 2009; Guillen et al., 2012; Datta et al., 2014; Tsakiroglou et al., 2015; Armstrong et al., 2016). A particular value in saturation does not necessarily imply that a unique disconnected structure (arrangement) of the nonwetting phase will settlein to comply with the externally imposed flow conditions, or (Ca, r), or − equivalently − by a set of two capillary numbers, (Ca_{o}, Ca_{w}), pertaining to oil and water (Tsakiroglou et al., 2015). Saturation provides only partial information about the process as − in essence − it represents an integral (volumetric) of an internal snapshot of the flow. The deficiency is not accounted on the integration per se (as long as − by definition − represents a macroscopic quantity) but on its inadequacy in describing (uniquely) the flow conditions − as it brings no definite input to the momentum balance. The issue is discussed in two recent papers (Valavanides et al., 2016, Section 2; Valavanides, 2018, submitted for review). Additional comments are provided in Section 4.9 ‘Discussion of results’ in the present work.
2 Pore Network
Implementation of the DeProF analysis and methodology presented in the previous paragraphs is directly related to the structure of the pore network and the associated modeling provisions pertaining to the prototype flows and, especially, to the disconnected fluidic elements. The following paragraphs describe the essential characteristics of the typical model pore network used in the present simulations.
2.1 Pore Network Geometry
The simulations are performed on the threedimensional (3D) cubic pore network models of the chamberandthroat type. This is a general type of model pore network; essentially it is a 3D version of the 2D network used in the laboratory studies of Avraam and Payatakes (1995, 1999).
The 3D cubic lattice pore network is of the “ballandstick” type, whereby spherical chambers, residing at the nodes of the cubic lattice, are interconnected with straight cylindrical throats, Figure 5. The network is constructed according to some typical distribution of the chamber and throat size (or size classes). The occurrence probabilities of size classes as well as the reduced dimensions of chambers and throats, of the 3D cubic lattice pore network used in the simulations, normalized to the lattice constant, are presented in Table 1. Each jikclass unit cell comprises two chambers of classes j and k with diameters D_{Cj} and D_{Ck}, interconnected with throati, a right circular cylinder of diameter D_{Ti}. The thick dashed line rectangle represents one of the three diagonal planes of the octahedron, while the dasheddotted lines are the two of the three octahedron diagonals (Fig. 5).
Occurrence probabilities, f_{j}, and respective reduced size diameters for the 5 classes of chambers, D_{C}_{j}, and throats, D_{T}_{j}, of the 3D pore network used in the DeProF simulations. All dimensions are normalized to the lattice constant, ℓ = 1221 μm.
A schematic representation of the network layout is shown in Figure 6(a). Chambers and throats are cited at the nodes and branches of the network. The mean coordination number is σ_{3D} = 6. For the particular network dimensions, the volume porosity is ϵ = 0.0464 and its absolute permeability is .
Figure 6 (a) A typical 3D cubic pore network. The lattice constant, ℓ, compared to the size of the throats and chambers, is shown expanded. The thick arrow indicates the direction of the macroscopic flow, (b) The pore network cubic lattice with one main diagonal parallel to the macroscopic flow direction. The two parallel equilateral triangles are perpendicular to the macroscopic flow, (c) Equilateral triangles align to form a checkerboard pattern and exactly fill any frontal area perpendicular to the macroscopic flow. 
There is no correlation between classes of chambers and throats therefore the network is isotropic on the macroscopic scale. Let ℓ denote the lattice constant. Each unitcell occupies an equilateral octahedron of pore network space with diagonal lengths equal to ℓ, edge lengths equal to and volume . The geometry of the pore space within every unitcell is defined by the shape of the chambers and throats (Fig. 5). For this particular type of cubic lattice networks, each unitcell comprises one eighth of each one of the two adjacent chambers that are interconnected with the throat.
An essential modeling consideration is that all fluidic elements move along the shortest of the pathways that are aligned to the macroscopic flow direction. That means that fluidic elements cannot occupy two parallel branches (throats) in the same lattice element (cube). For example, in Figure 6(a) a ganglion of size 6 is depicted as it is aligned to the macroscopic flow. The dotted line extending the ganglion ends outlines its trajectory. Any crosssectional plane, perpendicular to the macroscopic flow direction, would intersect any fluidic element or any fluid trajectory only once within a single throat or a single chamber.
The structure of this network is such that unit surface elements − equilateral triangles − completely fillup any crosssection perpendicular to the macroscopic flow (frontal area). Focusing on the macroscopic flow direction in Figure 6, in every unitlattice there exist two parallel equilateral triangles, oriented perpendicular to the main diagonal and to the macroscopic flow. These are depicted in Figure 6(b), one in plain light green color in front and the other, behind it, depicted with a shingled pattern. The triangles have side lengths equal to and surface area equal to . When the cubic lattice elements are stacked, these equilateral triangles stick sidebyside and exactly fill the frontal area. An impression of a 3×2×1 stacking of unitlattice elements (cubes) is depicted in Figure 6(c); cubes are aligned so that the main diagonal is parallel to the macroscopic flow direction and perpendicular to the figure. Five parallel layers, of stacked lattice surface elements, each comprising 1/3/4/3/1 equilateral triangles, are drawn one behind the other. Each layer exactly covers the respective frontal area. Along the main diagonal of the pore network, the zigzag corridors (or pathways), that are formed by concatenations of empty unitcells aligned to the macroscopic flow, are virtually confined within a normal equilateral twisting prism. The prism sides are equal to and its frontal area is equal to .
An important parameter associated with the ‘zigzag' motion of the ganglion is the tortuosity of its spine, χ_{G}. This is also the ratio of the length of its actual migration path over the effective migration length in the macroscopic flow direction. This parameter is essentially one of the p.m. geometrical and topologic characteristics and will be considered to be equal for all ganglion size classes. For the case under consideration, .
2.2 Types and number of unit cells in the Disconnected Oil Flow domain (DOF)
Ganglia belonging in the same sizeclass are considered to be indistinguishable. Nevertheless, the unit cells that they occupy, can be classified into four types, depending on whether these unit cells are conductive or blocked and on whether they host an o/w interface or not. The classification is used in accounting the surface area pertaining to the mobile /immobile o/w menisci. All droptraffic flow cells are identical (momentumwise) since they contain water and uniformly dispersed oil droplets.
Ganglion unit cells
The number of unit cells occupied by each iclass ganglia are given by (18) where N_{d} is the dimensionality of the network lattice. For the 3D lattice used in the simulations N_{d} = 3 (but for the 2D depiction of the lattice in Figure 3, N_{d} = 2)
Unit cells in which ganglion menisci exist are called gate unit cells. In the present work, we assume that all ganglia attain a ‘wormlike’ alignment − as much as the 3D zigzag pore conduits allow − to the macroscopic flow direction. Among all gate unit cells of a mobilized ganglion, there are only two gate unit cells in which rheons take place (one cell for the ‘xeron’ and one cell for the ‘hygron’ (Payatakes, 1982). Such cells are denoted with the index ‘E, G’ (endgate, ganglion), see Figure 7. The only potential endgate unit cells are those located at the foremost part (for the hygron) and at the aftmost part of the ganglion (for the xeron). All cells of any mobilized ganglion that are completely saturated with oil will be called ‘core cells', and will be denoted with the index ‘C,G’ (core, ganglion). The rest of the ganglion cells (which are neither endgate cells nor core cells) contain oil and water separated by virtually immobile menisci and are called sidegate unit cells. Ganglion migration induces a corresponding flow rate, , in the core cells and an equal flow rate in the endgate cells. This flow must be consistent with the mean velocity of the mass center of the iclass ganglion, , i.e., (19) where χ_{G} is the network tortuosity, defined in the previous subsection.
Figure 7 A closeup sketch view of Figure 1(b) indicating the different types of unit cells (u.c.) in the DOF domain. The classification of the ganglion unit cells (‘E,G’, ‘X,G’ and ‘C,G’) is according to their conductance state. 
The rest of the cells of a mobilized ganglion (sidegate unit cells) are considered as nonconducting, and they are denoted with the index ‘X,G’ (‘nonconducting, ganglion’). All of the cells of any iclass stranded oil ganglion are ‘X,G’ cells.
The fractional distributions of the various flow configurations of any iclass ganglion cells (mobilized as well as stranded ganglia taken into account), as indicated in Figure 7, are explicitly given by [ref. Equation (18) in Valavanides and Payatakes, 2001] as follows:
Core ganglion cells − conducting, (20) Nonconducting ganglion cells (note here, that this fraction takes into account all the cells of stranded ganglia), (21)End ganglion cells − conducting, through which hygrons take place, (22)Now, as expected, (23)
Drop Traffic Flow unit cells
The unit cells in the DTF domain contain water and oil droplets. Water may be connected or not depending on the formation of o/w menisci in contact with the pore walls. This is a modeling provision, associated − in general − to the wetting characteristics of the system, i.e. the combination of the dynamic contact angles, the size of the droplets and the pore wall microroughness. A DeProF model provision is that throats give birth to oil droplets having critical dimension (diameter) equal to the critical dimension (diameter) of the parent throat. Therefore, there are oil droplet size classes on a onetoone correspondence with pore/throat size classes. Nevertheless, as the population of oil droplets migrates downstream, they are mixed and there is no correlation between poresize and resident dropletsizes. The flow setup in the DTF domain unitcells is considered similar to oilinwater emulsion flow in capillary conduits (Cobos et al., 2009). The DTF unit cells are all considered to be indistinguishable − refer to Section 2.2 and Figure 7.
2.3 Mobilization of drops and ganglia
Oil droplets and ganglia have a tendency to get stranded, due to the capillary pressure difference that can be maintained across their o/w interfaces. This is attributed to the hysteresis between the advancing/receding contact angles that, in combination with the local pore geometry configuration, result in the hysteresis of the mean curvatures between the advancing and receding menisci. The mobilization probability of a ganglion or a droplet depends on the microscopic pressure difference − locally induced − along every pair of o/w interfaces and the net capillary pressure differences over the respective o/w interfaces (Payatakes, 1982).
In the DeProF model, stranding of ganglia and droplets is examined over particular stranding configurations, i.e. layouts within the unit cells they reside in, for which they get stranded. For these stranding configurations within the unit cells the capillary pressure differences become large and stable. Stability means that if any oil blob (drop or ganglion) is slightly perturbed, the new capillary pressure differences become larger and force the blob to return into its original position, keeping it stranded, until − eventually − another, larger pressure difference (over) perturbs the blob and destabilizes it; then it migrates downstream in consecutive episodic rheons, triggering Haines jumps and propagating pressure pulses (Berg et al., 2013; Rücker et al., 2015). For each type of blob and possible configuration, a respective mobilization probability is assigned to a solitary blob, i.e. a single blob in the pore network. Then, to account the crowding effect of the motion of neighboring blobs (whether drops or ganglia) induced by their episodic pore scale events, the locally induced macroscopic pressure gradient is perturbed by imposing a lognormal distribution around its mean value (Valavanides and Payatakes, 2001). This, on its turn, results in broadening the pressure interval from onset to complete mobilization of oil blobs (droplets and ganglia). Oil blobs as members of dense populations, have a higher mobilization probability than solitary ones.
The configuration for which a ganglion is least expected to get mobilized (or most probable to get stranded) is depicted in Figure 8a.
Figure 8 (a) A ganglion of size class 4 at a least expected to get mobilized configuration in the time interval between two rheons. (Emphasis is on the endgate unit cells), (b) The two configurations for which an oil droplet may become stranded. In both cases, the driving force for mobilization is from left to right. represents the idealized radii of curvature of the o/w interfaces for advancing and receding contact angles. 
A drop is least expected to get mobilized (most probable to get stranded) over two configurations, as outlined in Figure 8b: when invading a throat (denoted as Configuration I) or when displaced inside a throat (Configuration II). In reality, blocking of a unit cell by a drop may happen due to drop straining − when the critical dimension (volumetric diameter) of a drop is larger than the diameter of the throat, or when two small drops intercept at the throat entrance (Cobos et al., 2009). Both phenomena are observed in the laboratory study of Avraam and Payatakes (1995; videotapes).
DeProF model estimates of the mean mobilization probabilities for dense populations of droplets and ganglia of various sizes, as a function of the reduced macroscopic pressure gradient, x, are presented in the diagram of Figure 9. The diagram pertains to the mobilization of average size droplets and typical ganglion sizes of reduced volume equal to 1, 3, 5, 7, 9. The mean mobilization probability is defined as the weighted average of the mobilization probabilities of all configurations. (‘Weighted average’ stands for ‘multiplying the mobilization probability by the respective occurrence probability and summingup over all configurations’.) Note that ganglia are more difficult to mobilize than drops. The onset of mobilization takes place at values of x around x = 10. Note, also, that the mobilization probability of droplets cannot fall below a minimum value, even at the lowend of the pressure gradient domain. This is attributed to the DeProF modeling provision that throats give birth to spherical oil droplets having diameter equal to the parent throat diameter, therefore, as there is no correlation in the occurrence of droplet size and throat size in the network, it is always probable that a small droplet enters a large throat.
Figure 9 Mean mobilization probabilities for various types and sizes of oil blobs, as members of a dense population in terms of the reduced macroscopic pressure gradient, x. Continuous red line shows the mobilization probability of droplets, m^{D}, discontinuous black lines the mobilization probability of ganglia of various size classes, . 
3 Variables describing the interstitial flow structure
An inherent characteristic of immiscible twophase flow in porous media, is the disconnection of the nonwetting phase (oil). Fluid immiscibility, wetting, and pore network tortuosity, in combination with pore scale hydraulic instabilities and propagating pressure pulses (created by Haines jumps), trigger detachment, breakup and pinchoff of the nonwetting phase. These phenomena are counterbalanced − to a certain extent − by collision and eventual coalescence of nonwetting blobs (oil ganglia and oil droplets). The combined effect of these events leads to the collective behavior of the disconnectedoil flow; the latter can be described with population balances (Payatakes, 1982; Valavanides et al., 1998).
Oil disconnection may be accounted for by different variables depending on the specificity level of the required information. Ganglion and droplet size distribution, o/w surface area, coefficient of oil fragmentation, ganglion and droplet number density (overall population p.u.p.m.v.), FAV etc. are indicative variables − in decreasing order of specificity as to the transmission of information from the pore − to the macroscopic − scaleup.
The following paragraphs provide the basic methodology implemented in the DeProF model algorithm for estimating average ganglion sizes, surface area of o/w interfaces in the disconnectedoil flow domains and their partitioning over the GD and DTF domains. All variables are accounted after a physically admissible solution (i.e. conforming with the externally imposed flow conditions) has been detected. A bottomup process is implemented, e.g. the surface area of the o/w interfaces residing in the various types of unit cells − hosting parts of ganglia or droplets − is accounted, then averaged and summedup over the various flow configurations and/or the blob populations (ganglia and droplets) or the entire connected or disconnected flow domains.
3.1 Mean ganglion sizes
The expected values of the mean ganglion size are computed after a physically admissible solution, specifically, a physically admissible ganglion size distribution, , has been detected. Mean ganglion sizes are not simple averages but number density weighted averages.
In particular, the mean ganglion size pertaining to the entire ganglia population, , is estimated as the ratio of the entire oil ganglia mass over the entire ganglia population, as (24) where is the volume of oil comprising any ganglion of classi or, equivalently, the volume of iclass ganglia unit cells.
Similarly, the mean ganglion size pertaining to the mobilized ganglia population, , is estimated as the ratio of the entire oil mass of mobilized ganglia over the population of mobilized ganglia, as (25)
Consequently, the mean ganglion size pertaining to the population of stranded ganglia, , is estimated as (26)
3.2 Interstitial layout of the flow in terms of oil disconnectedness and o/w surface areas
The o/w surface areas are computed per unit volume of porous medium (p.u.v.p.m.), , to provide a rough indication of the structure of the entire flow, or per unit volume (p.u.v.) of the disconnected − oil flow domains, , e.g. p.u.v. of GCD, , or p.u.v. of DTF domain, . Still more information can be retrieved by accounting the specific (p.u.v.) o/w surface area pertaining to mobilized disconnectedoil blobs, , these being ganglia, , and/or droplets, , so that and , or pertaining to stranded oil blobs, . For simplicity, such estimates are not presented here.
Specific (p.u.v.p.m.) surface areas of o/winterface in the GCD region
Every mobilized ganglion migrates downstream through consecutive episodes of a xeron and a hygron (both called rheons). During a xeron, a receding o/winterface gets displaced across the frontmost gate unit cell hosting the head (bow) of the ganglion (Fig. 7). During a hygron, an advancing o/w interface gets displaced across the aftmost gate unit cell hosting the tail (stern) of the ganglion. It is only through those two gate unit cell (E,G) that the o/w interfaces move as the ganglion migrates downstream. The core unit cell (C,G) are completely saturated with oil and they host no o/w interface. On average, the mobile o/w surface area of any iclass ganglion is denoted as . This mobile o/w interface comprises exclusively the respective receding and advancing o/w interfaces in the two gate unit cells (E,G). The meanacrossthe gateu.c. surface of o/winterfaces, and , are computed for all jikclasses of unit cells (combinations of j and k chambers connected with ithroats) during the initial preprocessing steps (in the DeProF simulations) and considering the respective − advancing and receding − contact angles.
The rest of the o/winterfaces at the side u.c. (X,G cells) of any ganglion remain stranded within the throats of the side u.c. of the ganglion irrespective of its mobility (mobilized or stranded) (Fig. 7).
The following surfaces can be estimated:
Mean surface area of the mobile o/w interfaces, residing at the E,G unit cells of any mobilized ganglion, whereby the mean is considered over all E,G unit cell classes, (27)
Mean surface area of the stranded o/w interfaces in the X,G u.c. of an iclass ganglion, whereby the mean is considered over all X,G unit cell classes, (28)
Total mean surface are of the o/w interface of an iclass mobilized ganglion (29)
Total mean surface area of the o/w interface of an iclass stranded ganglion, whereby the mean is considered over all unit cell classes, (30)
In the previous expressions, , respectively represent the stranded o/w surface areas over all the classes of chambers (advancing interface) and throats (receding interface) considering the stranding configuration of any ganglion (Fig. 8).
One may now proceed with the computation of average macroscopic interfacial quantities, expressed per unit volume (p.u.v.) of the respective pore network cell domains, i.e. ganglion o/w interfaces in the Ganglion Cells Domain (GCD), , as well as for o/w interfaces in the DTF domain, A_{owMD}, A_{owSD}, A_{owD}. In all relevant expressions, the index notation consistently determines the physical quantities as follows: G – ganglia, D – Droplets, (ow)MMG/SMG – mobilized/stranded interface of mobilized ganglia, (ow)MD/SD – mobilized/stranded interface of droplets, (ow)G – total interface area of ganglia etc. In particular, the mobilized o/w surface area of mobilized ganglia p.u.v. of GCD is estimated by (31) the stranded o/w surface area of mobilized ganglia p.u.v. of GCD is (32)the total o/w surface area of mobilized ganglia p.u.v. of GCD is (33)
Similarly, the o/w surface area of stranded ganglia p.u.v. of GCD, is (34)
The stranded o/w surface area of ganglia (mobilized and stranded) p.u.v. of GCD, is (35)
The total o/w surface area of ganglia (mobilized and stranded) p.u.v. of GCD, is (36)
Specific (p.u.v.p.m.) surface areas of o/winterface in the DTF region
The average (mean) oil drop volume is (37)
The average number of oil drops p.u.v. of DTF domain is (38)
The o/w surface area of mobilized drops p.u.v. of DTF domain is (39) while the o/w surface area of stranded drops p.u.v. of DTF domain is (40)where and are the o/w surface areas of mobilized and stranded drops respectively, estimated as statistical averages over the various mobilization and stranding configurations. In particular, the mean o/w surface of mobile drops, , is computed as drops attain all possible configurations during their migration. Every nclass drop is screened and then its mean o/w surface area is computed as it crosses a jikclass unit cell. Similarly, the mean o/w surface area of stranded drops, , is computed as drops attain all possible configurations during their stranding. Every nclass drop is screened to compute its mean surface of the o/w surface area as it blocks a jikclass unit cell.
The total o/w surface area of (stranded and mobilized) drops p.u.v. of DTF domain can be estimated using Equations (39) and (40) as (41)
The total o/w surface area p.u.v. of DOF domain is (42)
The total o/w surface area p.u.v. of the whole pore network is (43)
Expected, ensemble average values of the reduced, specific o/w surface area in the DOF domains, for different values of the capillary number, the flow rate ratio and the viscosity ratio, are presented in Section 4.
Coefficient (fraction) of oil fragmentation
A measure of the extent of oil disconnection is provided by the coefficient of oil fragmentation, , defined as follows.
Consider the number density of the oil drops that could be formed if all of the oil in place had shuttered into oil droplets of minimum size, i.e. droplets with diameter equal to the minimum throat diameter. As already referenced in (2.2), throats give birth to spherical oil droplets having diameter equal to the parent throat diameter. Therefore, there are oil droplet size classes on a onetoone correspondence with thoat size classes. The minimum size of an oildroplet is, (44) whereas the corresponding o/w surface area is (45)
Now, consider the number density of the oil drops that could be formed if all of the oil in place had shuttered into oil droplets of minimum size, , (46)
The respective specific o/w surface area, , would comprise spherical surfaces pertaining to the mimimum size drops, of diameter equal to the diameter of the minimum throat, , therefore, (47) The coefficient (or, fraction) of oil fragmentation is defined as (48)
It is straightforward to derive its limiting values as , and based on its definition, the coefficient of oil fragmentation, f_{OF}, could also be called ‘interfacial saturation’. Expected, ensemble average values of are presented in Section 4.
3.3 Interstitial flow structure in terms of transport properties
To map the flow structure in terms of transport properties and not just the flow layout, one should account for dynamic variables that can efficiently describe the transport characteristics of the flow. In particular, ganglion and droplet velocity distribution, fractional flows of connected (CPF) and disconnectedoil flow (through GD and DTF), average size of mobilized and stranded ganglia and droplets, mobilized o/w surface area and its fractional decomposition over the GD and DTF domains, can describe the transport properties of disconnected flow − on different levels of specificity.
Interstitial fractional flows
The fractions of the oil flow rate in the CPF, GD and DTF domains are defined as, (49) where and are the reduced superficial velocities of oil in the GD and in the DTF domains respectively. Both share a common value with the reduced superficial velocity of oil in the disconnectedoil flow domain, , as the former comprise an inseparable mixture within DOF. They are given by (50)where is the reduced superficial velocity of oil in the (CPF), provided by Equation (6).
Expected, ensemble average values of the fractional flows through connected pathways, η_{o,C}, GD, η_{o,G}, and DTF, η_{o,D}, for different values of the capillary number, the flow rate ratio and the viscosity ratio, are presented in Section 4.
Superficial velocity of o/w surface in the GD domain
The superficial velocity of the o/w interface of the entire swarm of ganglia,, may be estimated as the total momentum of ganglia o/w surface divided by the total o/w surface area, (51) and then, reduced by the macroscopic velocity of oil, it becomes: (52)
Superficial velocity of o/w surface in the DTF domain
The reduced superficial velocity of the o/w surface in the DTF domain is readily estimated as (53)
It is worth noting the difference between the expressions in Equations (52) and (53). The summation in Equation (52) is to account the different momenta of the ganglion size classes migrating downstream at different velocities (to each ganglion size class corresponds a size class velocity), whereas in the indistinguishable DTF unit cells all of the drops that are mobilized migrate downstream at the uniform superficial velocity, .
Transport of o/w surface in the DOF domain
The reduced superficial velocity of the o/w interface in the DOF domain is estimated as (54)
Fractional o/w surface transport through GD and DTF
The fractional transport of the o/w interfaces through GD and DTF, and , may be estimated as the momentum of ganglia o/w surface divided by the total o/w surface area, (55)
Expected, ensemble average values of , and for different values of the capillary number, the flow rate ratio and the viscosity ratio, are presented in Section 4.
4 Simulations of steadystate 2ph flow in pore networks
A large number of simulations of steadystate twophase flows in the threedimensional (3D) chamberandthroat type pore network described in Section 2 were carried out implementing the algorithm build around the DeProF model. The values of the physicochemical properties of the pair of fluids (oil and water) examined in the simulations are: dynamic viscosity of water, , oil/water interfacial tension , advancing (A) and receding (R) oilwater/p.m. contact angles, and . Five typical systems were examined: favorable, κ = 0.33, 0.67, indifferent κ = 1.00 and unfavorable, κ = 1.50 and 3.00, viscosity ratio values were examined.
The simulations covered a rectangular domain logCa × logr, from (logCa, logr) = (− 8, − 2) to (logCa, logr) = (− 3, + 2). The domain was covered in 9 successive steps per order of capillary number magnitude, i.e. Ca = 1, ..., 9 × 10^{−j}, j = − 8, ..., − 4 (9 × 5 =45 divisions in the logCa range), and in 21 successive steps of logr = 0.2 in the logr range. Therefore, 45 × 21 = 945 simulations, pertaining to different flow conditions, were run for each of the 5 different values of the viscosity ratio. Each diagram in Figure 10 and in 12–19 contain the points that depict the average values of physical quantities (dependent variables) describing the interstitial structure of the flow for the corresponding set of flow conditions (Ca, r). Each diagram in these figures should normally contain 945 points. Nevertheless, this number is not reached as there are many flow conditions for which no solutions could be detected, i.e. steadystate twophase flow is not possible/sustainable. Because of space limitations, albeit 5 different viscosity ratio values were examined, only the results pertaining to κ = 0.33, κ = 1.00 and κ = 3.00 are presented. The flow structure maps of the other two systems (κ = 0.67 and 1.50) fall in between.
Figure 10 Expected, ensemble average values of the reduced pressure gradient, x, for different values of the capillary number, Ca, flow rate ratio, r, and viscosity ratio , Arrows indicate increasing Ca values: →10^{−8} ≤ Ca ≤ 10^{−4} and    > 4 × 10^{−7} ≤ Ca ≤ 10^{−4}. 
Figure 12 Expected, ensemble average values of the flow arrangement variables for different values of the capillary number, Ca, and viscosity ratio, κ. (a) Saturation, S_{w}; (b) volume fraction of connected pathway flow (CPF), β; (c) volume fraction of ganglion dynamics flow (GD) within the disconnected oil flow (DOF), ω. Arrows indicate increasing Ca values: →10^{−8} ≤ Ca ≤ 10^{−4} and    > 4 × 10^{−7} ≤ Ca ≤ 10^{−4}. 
Figure 13 Expected, ensemble average values of the reduced mechanical power dissipation, W, and its partition into the individual interstitial, constituent flows, for different values of the capillary number, Ca, flow rate ratio, r, and viscosity ratio, . Arrows indicate increasing Ca values: →10^{−8} ≤ Ca ≤ 10^{−4} and    > 4 × 10^{−7} ≤ Ca ≤ 10^{−4}. 
Figure 14 Expected, ensemble average values of the mean ganglion sizes for different values of the capillary number, Ca, and viscosity ratio, κ. a) Mean size of all ganglia (mobilized + stranded), ; b) mean size of mobilized ganglia, . Arrows indicate increasing Ca values: →10^{−8} ≤ Ca ≤ 10^{−4} and    > 4 × 10^{−7} ≤ Ca ≤ 10^{−4}. 
Figure 15 Expected, ensemble average values of the reduced o/w interface surface area for different values of the capillary number, Ca, flow rate ratio, r, and viscosity ratio, . (1st row) ganglion dynamics, A_{ow,G}; (2nd row) oil droplets, A_{ow,D}; (3rd row) Disconnected oil flow, A_{ow}; (4th row) overall. Arrows indicate increasing Ca values: →10^{−8} ≤ Ca ≤ 10^{−4} and    > 4 × 10^{−7} ≤ Ca ≤ 10^{−4}. 
Figure 16 Expected, ensemble average values of the coefficient (fraction) of oil fragmentation, f_{OF}, (or interfacial saturation) for different values of the capillary number, Ca, flow rate ratio, r, and viscosity ratio, . Arrows indicate increasing Ca values: →10^{−8} ≤ Ca ≤ 10^{−4} and    > 4 × 10^{−7} ≤ Ca ≤ 10^{−4}. 
Figure 17 Expected, ensemble average macroscopic fractional flows through connected pathways, η_{o,C}, ganglion dynamics, η_{o,G}, and drop traffic, η_{o,D}, flows. Arrows indicate increasing Ca values: →10^{−8} ≤ Ca ≤ 10^{−4} and    > 4 × 10^{−7} ≤ Ca ≤ 10^{−4}. 
Figure 18 Expected, ensemble average values of the o/w surface transport variables for different values of the capillary number, Ca, flow rate ratio, r, and viscosity ratio, . a) Reduced superficial velocity of o/w surface, U_{ow,DOF}, b) Fraction of o/w surface transport through ganglion dynamics , c) Fraction of o/w surface transport through drop traffic flow,, Arrows indicate increasing Ca values: →10^{−8} ≤ Ca ≤ 10^{−4} and    > 4 × 10^{−7} ≤ Ca ≤ 10^{−4}. 
Figure 19 Effect of the flow rate ratio on the energy efficiency of the process (oil output per kW of power dissipated by the process), f_{EU} = r/W, for different values of the capillary number, Ca, flow rate ratio, r, and viscosity ratio, . Arrows indicate increasing Ca values: →10^{−8} ≤ Ca ≤ 10^{−4} and    > 4 × 10^{−7} ≤ Ca ≤ 10^{−4}. 
The afore mentioned diagrams (Fig. 10, 12–19) share common characteristics. The abscissa axis measures the logarithm of the flow rate ratio. In every diagram, isocapillary points are outlined with curves. As the level of Ca increases in a stepwise manner (Ca = 1, ..., 9 × 10^{−j}, j = − 8, ..., − 4), the corresponding isocapillary curves shift from right to left and extend over narrower logr domains. The isocapillary curves for Ca = 10^{−j}, j = − 8, ..., − 4 are designed thicker for quick reference on changes in order of magnitude. In each diagram, arrows indicate the direction of arrangement of the isocapillary curves for increasing Ca values; a solid line arrow indicates changes over 10^{−8} ≤ Ca ≤ 10^{−4}, while a dashed line arrow is for changes over 4 × 10^{−7} ≤ Ca ≤ 10^{−4} (see point 2 below).
General remarks:

For all the systems examined, the higher Ca value for which the algorithm could detect a solution is Ca = 1.0 × 10^{−4}

For equal viscosity systems (κ = 1), the DeProF algorithm failed to detect any physically admissible solutions in the low Ca regime; specifically no solution could be detected for Ca < 4 × 10^{−7}. Presumably, viscosity disparity allows the process more degrees of freedom in adjusting itself within a richer ensemble of equiprobable physically admissible flow configurations (see discussion, 4.10).
The following Sections 4.1–4.8 provide an overall (basic) description of the furnished diagrams. Extended discussion of the particular trends in the flow arrangement/structure, as observed across flow conditions and system configurations, are provided in Section 4.9.
4.1 Reduced pressure gradient
Expected, ensemble average values of the reduced pressure gradient, x, Equation (5), for different values of the capillary number, Ca, and the flow rate ratio, r, are presented in the diagrams of Figure 10 for different viscosity ratio systems, κ.
Observations of single phase flows within pore networks confirm that the macroscale pressure gradient scales linearly with the superficial fluid velocity governed by Darcy's law. This seems to be a quite trustworthy modeling consideration also in the case of twophase flow but only when very high superficial velocities are examined, whereby the capillary forces are negligible. However, at moderate/low velocities, when the capillary forces are comparable to the viscous forces, the macroscopic pressure gradient in the steadystate does not scale linearly with the flow rate.
The DeProF model predictions in Figure 10 are consistent with recent experimental studies on steadystate twophase flows in glass beads in HeleShaw flow experiments (Tallakstad et al., 2009), in virtual flows in 3D pore networks (Sinha and Hansen, 2012), in experiments in glass bead columns (Sinha et al., 2017), in 2D glassetched pore network models and in sandpack column experiments (Tsakiroglou et al., 2015), observing that the nonlinear relation between pressure gradient and the flow rate obeys generic power laws with different exponent values. Differences between the values of the scaling exponents in these studies attributed to differences associated to (a) dimensionality of the pertinent variables (b) measurements pertaining to different flow conditions, (c) dimensionality of the o/w/pm system properties. The DeProF model predictions bridge these discrepancies and show the potential for deriving a universal power law between appropriate dimensionless variables.
4.2 Typical, physically admissible flow configurations
Indicative physically admissible flow configurations described by the values of the FAV, , and the associated values of the ganglion size distribution parameters, , for different flow conditions, are presented in the diagrams of Figure 11. ‘Moderate flow arrangements’ are characterized by physically admissible flow arrangement values close to the ensemble average, (S_{w}, β, ω) and are presented in diagrams (a1)–(c1). ‘Extreme flow arrangements’ are characterized by physically admissible flow arrangement values pertaining to the extreme corners of the ensemble envelope and are presented in diagrams (a2)(c2).
Figure 11 Diagrams of the ganglion size distributions, pertaining to moderate physically admissible flows for different macroscopic flow conditions, (a1), (b1), (c1), and to extreme physically admissible flows for the same macroscopic flow conditions (a2), (b2), (c2). Tabulated data indicate the values of flow arrangement variables, and of corresponding ganglion size distributions parameters, . The values of flow conditions, (Ca, r), and average flow arrangement variables,(S_{w}, β, ω), are provided in the upper two rows of each tablet. 
4.3 Flow arrangement variables (steadystate average)
Expected values of the FAV, i.e. the saturation, S_{w}, the volume fraction of CPF, β, the volume fraction of GD within the DOF, ω, for different values of the capillary number, Ca, and the flow rate ratio, r, for different viscosity ratio systems, κ, are presented in the diagrams of Figure 12. As the viscosity ratio is increased, the onset and buildup of connected oil pathways (β) are observable from lower Ca values.
4.4 Reduced mechanical power dissipation
The reduced rate of mechanical energy dissipation within the process, W, is defined as (Valavanides et al., 2016)
(56)
where is the mechanical power dissipation for the equivalent onephase flow of water with superficial velocity, . The mechanical power dissipated in Connected Pathway Flow (CPF), Ganglion
Dynamics (GD) and DropTraffic Flow (DTF) are given respectively by
(57)
Expected values of the reduced mechanical power dissipation, W, and its partition into the individual interstitial, constituent flows, for different values of the capillary number, Ca, and the flow rate ratio, r, for different viscosity ratio systems, κ_{,} are presented in the diagrams of Figure 13.
4.5 Mean ganglion size (overall and mobilized)
Expected, values of the mean ganglion sizes, averaging Equations (24) and (25) over the ensemble of physically admissible solutions, for different values of the capillary number, Ca, and viscosity ratio, κ, are presented in the diagrams of Figure 14. The mean size pertaining to the entire ganglia population (mobilized and stranded), , is presented in (a), whereas the mean size of the mobilized ganglia,, in (b).
4.6 Interfacial areas and oil disconnectedness
Expected, ensemble average values of the reduced o/w interface surface area for different values of the capillary number, Ca, flow rate ratio, r, and viscosity ratio, are depicted in the diagrams of Figure 15. The specific o/w surface areas attributed to ganglia, A_{ow,G}, and oil droplets,A_{ow,D}, are provided on the 1st and 2nd row diagrams respectively.
As expected, the o/w surface area density in the DTF domain is much larger than that in the GD domain − in the cases considered here, 10fold higher. The diagrams on the 3rd row indicate the expected values for the o/w surface area density in the DOF domain, A_{ow,DOF}, as a weighted sum of A_{ow,G} and A_{ow,D}, see equation (42). The 4th row diagrams provide the expected values of the o/w surface density over the whole pore network volume A_{ow}, see equation (43). These diagrams provide a rough indication of the structure of the flow and can also decompose into other constituents, e.g. o/w surface area densities pertaining to mobilized disconnected oil blobs, , whether these are ganglia, , or droplets, , so that and . The significant drop in the overall o/w interface area density, A_{ow}, observed in Figure 15d at intense twophase flow conditions (high Ca, high r values), is due to the increase of the volume fractions of CPF, see Figure 12b.
The extent of oil disconnectedness can be accounted for by the coefficient of oil fragmentation, f_{OF}, Equation (48). Expected values of f_{OF}, for different values of the capillary number, Ca, flow rate ratio, r, and viscosity ratio, κ, are presented in the diagrams of Figure 16. In particular, f_{OF} can be used to assess the unattained margin of o/w surface area density. It can be considered as a reduced A_{ow}, normalized by the maximum value of A_{ow} that would be attained if all of the oil would shutter in droplets of the smallest possible size. In that context, f_{OF} could also be called ‘interfacial saturation’. As the flow rate ratio, r, increases, f_{OF} drops significantly, indicating that a part of DOF transforms (reorganizes) into connectedoil pathways (no o/w interface). This drop becomes more intense by shifting into higher isocapillary curves, as observed by the steeper drop from higher values in f_{OF} as flow conditions shift to higher isocapillary curves. The explanation is that as Ca increases, i.e. superficial velocity increases, the flow of DOF shifts mostly into DTF, while large ganglia dynamics turn to small ganglia dynamics. Interestingly so, the same transition is observed in S_{w}, Figure 12b.
4.7 Dynamic variables describing the interstitial transport properties of the process
Ensemble average values of the fractional flow of oil through CPF, η_{o,C}, GD, η_{o,G}, and DTF, η_{o,D}, Equations (49), are presented in the diagrams of Figure 17, for different values of the capillary number, the flow rate ratio and the viscosity ratio.
In all cases examined, the contribution of disconnectedoil flow to the total flow is substantial. Moreover, for some combinations of flow conditions (Ca, r) and viscosity ratio, DOF prevails. In DisconnectedOil Flow (DOF), the contribution of DTF is generally more intense than that through GD. At relatively lowCa/lowr values the fractional flow of oil through connected pathways is practically zero and most of the oil flow is disconnected − mainly in the form of (mobile) droplets and (less mobile) ganglia. Only when r is increased does CPF takesover a substantial fraction of the oil flow over droplets and ganglia. Then, with increasing Ca, this takeover shifts to lower r values and most of the oil flows through connected pathways.
Ensemble average values of the transport of o/w interfaces in DOF, U_{ow,DOF}, as well as the relative contributions in o/w transport through DTF, ξ_{ow,D}, and GD, ξ_{ow,G}, Equations (54) and (55), are presented in the diagrams of Figure 18, for different values of the capillary number, the flow rate ratio and the viscosity ratio.
In all cases, the contribution of transport of o/w interface through DTF, prevails over GD; in particular, across the lowend domain of flow rate ratio, r, the o/w interface transport takes place almost exclusively through DTF. This is attributed to the large volumetric density of o/w interfaces in DTF, Figure 15, in combination to the fractional flow of disconnected oil through DTF, Figure 17, across the entire domain of flow conditions.
4.8 Energy efficiency of steadystate twophase flow in porous media processes
The efficiency of any steadystate twophase flow in porous media process, with respect to the maximization of the oil transport per kW of mechanical power supplied to it (or ‘oil output per unit energy cost’), may be assessed by the values of the energy utilization coefficient, f_{EU}, a macroscopic quantity originally defined by Valavanides and Payatakes (2003) [see also Valavanides, 2012; Valavanides et al., 2016, Eqs. (6)–(8)] as (58)
Simulations implementing the DeProF model, suggest that a set of flow conditions exist for which the sought process attains locally maximum energy efficiency, see Figure 19. Such conditions, define a locus of optimum operating conditions or critical flow conditions, r*(Ca), for each oilwaterp.m. system (Valavanides, 2014 and Valavanides, 2018). The locally maximum value in f_{EU,ij} (Ca_{i}, r_{j}) detected in the simulations per fixed Ca_{i}, corresponds to a particular flow rate ratio value, .The sets detected for each system are depicted in the diagrams of Figure 20. The corresponding markers show a well defined trend: above a certain value in Ca the flow conditions for maximum efficiency become insensitive to Ca, indicating a change in the character of the flow, from capillarity − to viscositydominated. This particular Ca value − segregating the two types of flow − depends on the system properties. Note that in the diagrams of Figures 19 and 20, the horizontal, dashdotted line in red indicates the theoretically predicted values of the (global) maximum attainable energy efficiency (the o/w/pm system's ‘efficiency ceiling’), , while the vertical, dashed line in red indicates the corresponding flow rate ratio, . These conditions are theoretically met for ‘infinite’ Ca, in which case, o/w interfaces have no effect at all to the viscosity dominated flow.
Figure 20 The set of critical flow conditions, , detected by the DeProF simulations, for which the energy efficiency takes maximum values (ref. Fig. 19). 
Detecting and setting critical flow conditions r * (Ca) (of locally maximum energy efficiency) in a real process could potentially provide large marginal benefits in industrial applications, such as enhanced oil recovery (Taber et al., 1997, 1 and 2; Kjarstad and Johnsson, 2009; Alvarado and Manrique, 2010). The existence of ‘optimum conditions’ for oil transport in steadystate twophase flow in pore networks is a consequence of the remarkable internal adaptability of the flow to externally imposed flow conditions (Ca, r), see diagrams in Figure 13 pertaining to the internal partitioning of the total power dissipation into the constituent interstitial flows, Equations (56) and (57), and its inherent characteristic in tradingoff between CPF, GD and DTF and selfadjusting the connected versus disconnected movingoil balance.
DeProF model predictions with regards to the existence of a locus of optimum efficiency conditions have been recently substantiated by a post examination of relative permeability data, contained within a variety of published relative permeability diagrams, pertaining to steadystate twophase flows in porous media179 diagrams in total − for a variety of systems and flow conditions (Valavanides et al., 2016). The existence of critical flow (optimum operating) conditions is therefore an inherent characteristic of steadystate twophase flows in real porous media, remaining in latency until indicated by DeProF model predictions (Valavanides and Payatakes, 2003). Recently, in order to develop a theoretical framework based on statistical thermodynamics that would justify the existence of a locus optimum operating conditions, Valavanides and Daras (2016) proposed a model for accounting the number of the process microstates.
4.9 Discussion of results
The diagrams presented in Figures 10–19 indicate that distinct mutations can be observed. Some mutations are obvious and measurable at the macroscopic scale, others remain in latency within the interstitial flow configurations and, in practice, they are very difficult to measure. Nevertheless, the DeProF model has the capability to predict what could be revealed in a systematic laboratory study, would the necessary equipment and resources be available.
The behavior of the steadystate twophase flow structure for different flow conditions may be revealed by a comparative observation of the diagrams in Figures 10, 12–19, as follows.
Low Ca − The reduced pressure gradient diagrams of Figure 10 indicate that at the low Ca domain, the flow rate ratio, r, has no effect on the reduced pressure gradient, x. The later attains practically constant x values, indicating a saturatedflow type response of the process to flow conditions: At low Ca flow conditions, disconnected oil remains stranded and blocks part of the p.m. cross section that would otherwise (if mobilized) be available to the flow. This raises the reduced pressure gradient by a factor 100× to even 1000× as Ca is reduced to very low values (4×10^{−7} for medium to 10^{−4} for very low Ca). Getting deeper in the interstitial flow configuration, a mutation can be revealed as the flow rate ratio is increased from relatively low to higher values. Specifically,
Low Ca − low r − At the low end of the flow rate ratio domain, the prevailing flow pattern is (DOF) as practically no (CPF), β → 0, can be detected in the diagrams of row (b), Figure 12. Delving deeper in DOF, DTF occupies a significantly larger − over GD − fraction of the available p.m. crosssection, as is indicated by the rather low ω and high (1 − ω) values predicted in the diagrams of Figure 12c. Still for relatively low r values, DTF is also the prevailing flow pattern as is indicated in the oil fractional flows diagrams, Figure 17, where CPF is practically nonexisting [diagrams in row (a)], and oil flow through DTF [row (c)] prevails over oil flow through GD [row (b)]. What are the causes of this behavior? Small ganglia are difficult to mobilize when compared to larger ganglia and droplets, Figure 9. Therefore it is these small ganglia that get stranded and block a great part of the available p.m. crosssection, whereas it is only with the larger, mobilized ganglia that oil transport takes place through GD. This is indicated in the diagrams of Figure 18, whereby the mean size of mobilized ganglia is [diagrams in row (b)], large enough, if compared to the mean size of all ganglia, [diagrams in row (a)]. These statements conform absolutely with the predicted values of the o/w surface area densities, Figure 15, of the coefficient of oil fragmentation, Figure 16, and most importantly, the obvious prevalence of DTF − over GD − as the main pattern of o/w interface transport, as shown in Figure 18, diagrams in row (b) against those in row (c).
Low Ca − high r − With further increase in the flow rate ratio (still in the low Ca domain), connected pathways form, while the DOF reconfigures itself to a prevailing DTF mode. Still the prevailing flow pattern is expected to be Large GD as indicated by the large average size of mobile ganglia, Figure 14b, Figure 15 and the corresponding relatively small values in oil fragmentation, Figure 16.
High Ca − At conditions of very high capillary number the reduced pressure gradient attains values very close to 1. It behaves like its equivalent saturated flow. This is an indication that the capillarity effects are negligible when compared to the effects of bulk viscosities. This does not mean there are no o/w interfaces; as it will be described in the following, the flow is rich in o/w interfaces. Nevertheless their contribution on flow impedance is not significant.
High Ca − transition from low to high r– The flow pattern is transitional across the flow rate ratio domain eligible for twophase flow (−2 < logr < 0), from pure DOF, (1 − β) → 1, and borderline CPF β → 0 for −2 < logr < − 1, to mixed DOF and CPF for −1 < logr < 0, see Figure 12b. Delving deeper in DOF, DTF occupies a predominant fraction of the available p.m. crosssection − over GD, as is indicated by the low ω and high (1 − ω) values on the diagrams of Figure 12c across the whole range of flow rate ratio values (−2 < logr < 0). Albeit CPF does not occupy a significant volume fraction (topping up to 0.22 for κ = 0.33), its contribution in the oil transport is dramatic as it suddenly stepsup and takes over more than 60% of the oil flow rate as the flow rate ratio is increased, Figure 17a. This mutation takes place at logr ≈ − 1.5 for κ = 0.33. With further increase in r the fractional flow of oil through CPF climbsup towards 100%, until twophase flow is no more sustainable. The transition just described shifts towards lower r values for increasing values of the viscosity ratio. Focusing on the partition of oil fractional flow through GD and DTF, no significant change can be observed, Figure 17b and c. What are the causes of this behavior? Small ganglia are difficult (less probable) to mobilize when compared to larger ganglia and droplets, Figure 9. The flow structure is rich in small ganglia , see Figure 14a. The mean size of mobilized ganglia is also small , Figure 14b, because of the large pressure gradients. Nevertheless, the low volumetric occupation of GD in combination with the low migration velocity of small ganglia (as the driving force − the externally induced net pressure drop − is small) result in small fractional flow through GD, as attested in Figure 17b. The fractional flow of oil through DTF (the GD complementary constituent in DOF) is significant; nevertheless, it drops dramatically as CPF takes over, as attested in the diagrams in Figure 17c and a respectively. These statements conform with the predicted values of the o/w surface area densities, Figure 15, of the coefficient of oil fragmentation, Figure 16, and the obvious prevalence of DTF as the main transporter of o/w interface over GD, as shown in Figure 18, diagrams in row (b) against those in row (c).
Volume fractions against fractional flows
The previous analysis leads to the issue of whether using volume fractions in describing the flow, is more specific than using fractional flows. If we focus on the description of connected oil flow (the CPF domain), we may observe significant discrepancies between volume fractions and fractional flows − already identified in the previous paragraphs. In general, for the same flow conditions and viscosity ratios, the predicted values of the fractional flow rates, η_{o,C}, are larger than the corresponding values of the volume fractions, β. This is expected if we consider that the absence of o/w interfaces in the CPF domain does not impede the flow (as the the bulk viscosity does). Let us compare the diagrams of the volume fraction of connectedoil pathways, β, Figure 12b, with those of the fractional flow of oil through CPF, η_{o,C}, Figure 17b. Indicative values are tabulated for typical flow conditions and viscosity ratios, see Table 2. The discrepancies are pronounced in the low Ca / low r regimes and they tend to fadeaway with increasing Ca and/or r. The same behavior is observed for changes in the viscosity ratio. As a general rule of thumb, the discrepancies between volume fraction and fractional flow become less pronounced with increasing flow intensity, i.e. high Ca and/or, high r and/or increasing viscosity ratio.
Volume fraction of the connected pathway flow (CPF) domain, β, against fractional flow of oil through it, η_{o,C}, for various flow conditions, (Ca, r), and viscosity ratios κ.
Albeit there is a significant and observable mutation in the interstitial structure of the flow, overall, the macroscopic flow is stable. Variables that describe the interstitial structure have a different behavior than variables describing macroscopic, averaged physical quantities; the former are obviously more sensitive to changes in flow conditions that the latter. Compare, for example, the diagrams of the reduced pressure gradient, x, in Figure 10 and the energy efficiency index, f_{EU}, in Figure 19, pertaining to the entire flow, against the diagrams in Figures 15, 17 and 18, pertaining to interstitial flow characteristics of the two subdomains of DOF, GD and DTF. The behavior indicated in some of the diagrams of Figures 12, 14 and 16, pertaining to DOF, fall in between the behavior observed in the two previous groups.
Effect of viscosity disparity
In all cases examined, there is an obvious, universal trend in the dependence of the variables on logr: the set of isocapillary curves shift and expand with viscosity disparity. This trend can be observed by comparing any triplet of diagrams on the same row in Figures 12–19. The effect of viscosity ratio has been studied by Vizika et al. (1994).
The present simulations show that, if the twofluids have equal viscosities, twophase flow confines itself within a narrow domain of flow conditions, pertaining to relatively high Ca values (compare middle column diagrams with left and right columns, Figure 10; also Figures 12–19). This behavior is attributed to an intrinsic, extra degree of freedom provided to the flow by viscosity disparity. In this context, it is interesting to examine the two cases. If the two fluids have equal viscosities the flow is regulated exclusively by interfacial phenomena, i.e. interfacial tension, wetting and advancing /receding contact angle hysteresis. The latter prevail in the lowCa regime, immobilizing a large part of the disconnected oil, so that a sub domain of − relatively high − oil/water flow rate ratio values are physically unattainable. In the highCa regime, mobilization is more probable, therefore twophase flow is possible for lowr values. In that case (κ = 1), the ensemble of physically admissible flow configurations will comprise interstitial flow arrangements that would be regulated by interfacial phenomena. If the two fluids have different viscosities, the flow is sustained over a broader Ca, r domain. This is explained as follows: Viscosity disparity has an asymmetric effect; the bulk viscosity stresses induced by strain rates can only be balanced by increased pressure gradient that, in turn, expedites mobilization of the nonwetting disconnected fluidic elements. It should be stressed here that − according to predictions − the outcome will depend on the level of disparity but not on the type of disparity, i.e. whether favorable or unfavorable viscosity ratios. The analysis is supported by indicative values tabulated for typical flow conditions and viscosity ratios, see Table 3.
Reduced pressure gradient values, x, for different flow conditions, (Ca, r), and viscosity ratios, κ.
The systematic behavior observed across similar diagrams for different values of the viscosity ratio, especially those diagrams pertaining to macroscopic quantities, suggest that these (diagrams) may eventually collapse into a universal diagram if the capillary number is appropriately reduced by the viscosity ratio. This is discussed in a paper by Valavanides (2018).
Conclusion
Extensive simulations implementing the DeProF mechanistic model were run to deliver maps of physical quantities that describe the interstitial flow structure of twophase flows in a typical 3D model pore network of the chamberandthroat type, over a broad domain of capillary number, flow rate ratio and viscosity ratio. The scope was to understand the mechanisms underlying the significant mutations observed in twophase flows in porous media in mixtures of connected and disconnected flows (connectedoil pathway flow, dynamics of small and large ganglia, drop traffic flow).
Systematic mutations have been revealed within the flow structure with changing flow conditions. Some of them surfaceup on the macroscopic scale and can be measured e.g. the reduced pressure gradient. Other remain in latency within the interstitial flow structure e.g. the volume fractions or fractional flows of oil through connecteddisconnected flows. Deeper within the disconnectedoil flow, strong/intense mutations between Ganglion Dynamics and DropTraffic Flow have also been identified.
Mutations shift and/or become pronounced with viscosity disparity. The mutations are more evident over variables describing the interstitial transport properties of process than variables describing volume fractions. Τhis characteristic behavior is attributed to the interstitial balance between capillarity and bulk viscosity. Flow of the nonwetting phase in disconnected form is significant and for certain flow conditions the dominant flow mode.
In general, sustainability of twophase flow is in favor of disparities in the flow conditions (at relatively large values of the flow rate ratio) as well in the system configuration, e.g. the viscosity ratios away from unity (κ ≠ 1, 0), increased contact angle hysteresis (this last statement is not examined in the present work but can be inferred with confidence). This is attributed to the effect of the interstitial balance between capillarity and bulk viscosity over asymmetrical effects induced by disparity in physicochemical properties of the two fluids.
Albeit significant mutation is observed within the interstitial structure of the flow, yet, overall, the macroscopic flow is stable. Variables that describe the interstitial structure are obviously more sensitive to changes in flow conditions than variables describing macroscopic, averaged physical quantities.
A general phenomenological behavior, the nonlinearity associated with the reconfiguration of the interstitial flow structure (mutation), is evident in the disconnected nonwetting flow regimes (GD and DTF). This can be attributed to the contribution of the o/w interfaces and the associated capillary effects that, together with viscosity disparity, regulate (in a sense, they litigate) the partition of the entire flow into connected and disconnected nonwetting phase flow.
Concluding, the present work indicates the potential of the DeProF truetomechanism theoretical framework in revealing the effects of underlying flow mechanisms in immiscible steadystate twophase flow in porous media processes, justify the phenomenology and resolve the conflicting observations on whether disconnected oil moves or not and at what degree, in a systematic and consistent approach. Even its findings are based on a mechanistic model (theory), it can provide input for designing a smart, systematic laboratory study (experiment) that would eventually validate its predictions and expose more underlying mechanisms.
Acknowledgment
This research work has been cofunded by the European Union (European Social Fund) and Greek national resources under the framework of the “Archimedes III: Funding of Research Groups in TEI of Athens” project of the “Education and Lifelong Learning” Operational Program (MlS379389).
Dr. Olga Vizika, IFP Energies Nouvelles, is greatly acknowledged for discussions on the subject. Mrs. Evangelia Ponta, TEI Athens, is greatly acknowledged for assisting with the simulations. Valuable remarks and constructive comments from anonymous referees are also gratefully acknowledged.
References
 AlGharbi M.S., Blunt M.J. (2005) Dynamic network modeling of twophase drainage in porous media, Phys. − Rev. E 71, 016308, doi: 10.1103/PhysRevE.71.016308. [CrossRef] [Google Scholar]
 Alvarado V., Manrique E. (2010) Enhanced oil recovery: an update review, Energies 3, 1529–1575, doi: 10.3390/en3091529. [CrossRef] [Google Scholar]
 Armstrong R.T., McClure J.E., Berrill M.A., Rücker M., Schlüter S., Berg S. (2016) Beyond Darcy's law: The role of phase topology and ganglion dynamics for twofluid flow, Phys. − Rev. E 94, 043113, doi: 10.1103/PhysRevE.94.043113. [CrossRef] [Google Scholar]
 Avraam D.G., Kolonis G.B., Roumeliotis T.C., Constantinides G.N., Payatakes A.C. (1994) Steadystate twophase flow through planar and nonplanar model porous media, Transp. Porous Media. 16, 75–101. [CrossRef] [Google Scholar]
 Avraam D.G., Payatakes A.C. (1995) Flow regimes and relative permeabilities during steadystate twophase flow in porous media, J. Fluid Mech. 293, 207–236. [CrossRef] [MathSciNet] [PubMed] [Google Scholar]
 Αvraam D.G., Payatakes A.C. (1999) Flow mechanisms, relative permeabilities and coupling effects in steadystate twophase flow in porous media case of strong wettability, Ind. Eng. Chem. Research 38, 778–786. [CrossRef] [Google Scholar]
 Bazylak A., Berejnov V., Markicevic B., Sinton D., Djilali N. (2008) Numerical and microfluidic pore networks: Towards designs for directed water transport in GDLs, Electrochimica Acta 53, 7630–7637, doi:10.1016/j.electacta.2008.03.078. [CrossRef] [Google Scholar]
 Berg S., Holger O., Klapp S.A., Schwing A., Neiteler R., Niels Brussee N., Makurat A., Leu L., Enzmann F., Schwarz JO., Kersten M., Irvine S., Stampanoni M. (2013) Realtime 3D imaging of Haines jumps in porous media flow, Proc. Nat. Academy Sci. 110, 10, 3755–3759, www.pnas.org/cgi/doi/10.1073/pnas.1221373110. [CrossRef] [Google Scholar]
 Cobos S., Carvalho M.S., Alvarado V. (2009) Flow of oil–water emulsions through a constricted capillary, Int. J. Multiph. Flow 35, 507–515, doi: 10.1016/j.ijmultiphaseflow.2009.02.018. [CrossRef] [Google Scholar]
 Constantinides G.N., Payatakes A. C. (1996) Network simulation of steadystate twophase flow in consolidated porous media, AIChE J. 42, 369–382, doi: 10.1002/aic.690420207. [CrossRef] [PubMed] [Google Scholar]
 Constantinides G.N., Payatakes A.C. (2000) Effects of precursorwetting films in immiscible displacement through porous media, Transp. Porous Media 38, 291–317, doi: 10.1023/A:1006557114996. [CrossRef] [Google Scholar]
 Datta S.S., Ramakrishnan T.S., Weitz D.A. (2014) Mobilization of a trapped nonwetting fluid from a threedimensional porous medium, Phys. Fluids 26, doi: 10.1063/1.4866641. [Google Scholar]
 Fusseis F., Xiao X., Schrank C., De Carlo F. (2014) A brief guide to synchrotron radiationbased microtomography in (structural) geology and rock mechanics, J. Struct. Geol. 65, 1–16, doi:10.1016/j.jsg.2014.02.005. [CrossRef] [Google Scholar]
 Georgiadis A., Berg S., Makurat A., Maitland G., Ott H. (2013) Porescale microcomputedtomography imaging: Nonwettingphase clustersize distribution during drainage and imbibitions, Phys. − Rev. E 88, 033002, doi: 10.1103/PhysRevE.88.033002. [CrossRef] [Google Scholar]
 Ghassemi A., Pak A. (2011) Numerical study of factors influencing relative permeabilities of two immiscible fluids flowing through porous media using lattice Botzmann method, J. Pet. Sci. Eng. 77, 135–145, doi: 10.1016/j.petrol.2011.02.007. [CrossRef] [Google Scholar]
 Guillen VR, Romero MI, Marcio da Silveira Carvalho MS, Alvarado V. (2012) Capillarydriven mobility control in macro emulsion flow in porous media, Int. J. Multiph. Flow 43, 62–65, doi: 10.1016/j.ijmultiphaseflow.2012.03.001. [CrossRef] [Google Scholar]
 Kirkpatrick S. (1973) Percolation and conduction, Rev. Mod. Phys. 45, 4, 574–588, doi: 10.1103/RevModPhys.45.574. [CrossRef] [Google Scholar]
 Kjarstad J., Johnsson F. (2009) Resources and future supply of oil, Energy Policy 37, 441–464, doi: 10.1016/j.enpol.2008.09.056. [CrossRef] [Google Scholar]
 Knudsen H.A., Aker E., Hansen A. (2002) Bulk flow regimes and fractional flow in 2d porous media by numerical simulations, Transp. Porous Media 47, 99–121, doi: 10.1023/A:1015039503551. [CrossRef] [Google Scholar]
 Knudsen H.A., Hansen A. (2002) Relation between pressure and fractional flow in twophase flow in porous media, Phys. Rev. E 65, 0563101–05631010, doi: 10.1103/PhysRevE.65.056310. [CrossRef] [Google Scholar]
 Krummel A.T., Datta S.S., Münster S., Weitz D.A. (2013) Visualizing multiphase flow and trapped fluid configurations in a model threedimensional porous medium, AIChE J. 59, 1022–1029, doi: 10.1002/aic.14005. [CrossRef] [Google Scholar]
 Nguyen V.H., Sheppard A.P., Knackstedt M. A., Pinczewski W. (2006) The effect of displacement rate on imbibition relative permeability and residual saturation, J. Petrol. Sci. Eng. 52, 54–70, doi: 10.1016/j.petrol.2006.03.020 [CrossRef] [Google Scholar]
 Oughanem R., Youssef S., Bauer D., Peysson Y., Maire E., Vizika O. (2015) A multiscale investigation of pore structure impact on the mobilization of trapped oil by surfactant injection, Transp. Porous Media 109, 673–692. [CrossRef] [Google Scholar]
 Pan C., Hilpert M., Miller C.T. (2004) LatticeBoltzmann simulation of twophase flow in porous media, Water Resour. Res. 40, W01501. [Google Scholar]
 Payatakes A.C. (1982) Dynamics of oil ganglia during immiscible displacement in waterwet porous media, Ann. Rev. Fluid Mech. 14, 365–393. [CrossRef] [Google Scholar]
 Ramstad T., Idowu N., Nardi. C., Øren PE. (2012) Relative permeability calculations from twophase flow simulations directly on digital images of porous rocks, Transp Porous Med 94, 487–504. [CrossRef] [Google Scholar]
 Rücker M., Berg S., Armstrong R.T., Georgiadis A., Ott H., Schwing A., Neiteler R., Brussee N., Makurat A., Leu L., Wolf M., Khan F., Enzmann F., Kersten M. (2015) From connected pathway flow to ganglion dynamics, Geophys. Res. Lett. 42, 3888–3894. [CrossRef] [Google Scholar]
 Ryoo S., Rahmani A.R., Yoon K.Y., Prodanovic M., Kotsmar C., Milner T.E., Johnston K.P., Bryant S.L., Huh C. (2010) Theoretical and experimental investigation of the motion of multiphase fluids containing paramagnetic nanoparticles in porous media, Paper SPE134879 presented at the SPE Annual Technical Conference and Exhibition, 19–22 September, Florence, Italy, 20p. [Google Scholar]
 Sahloul N.A., Ioannidis M.A., Chatzis I. (2002) Dissolution of residual nonaqueous phase liquids in porous media: porescale mechanisms and mass transfer rates, Adv. Water Resour. 25, 1, 33–49, doi: 10.1016/S03091708(01)000252. [CrossRef] [Google Scholar]
 SerresPiole C., Preud'homme H., MoradiTehrani N., Allanic C., Jullia H. Lobinski R. (2012) Water tracers in oilfield applications: Guidelines, J. Pet. Sci. Eng. 9899, 22–39, doi: 10.1016/j.petrol.2012.08.009. [CrossRef] [Google Scholar]
 Sinha S., Hansen A. (2012) Effective rheology of immiscible twophase flow in porous media, Europhysics Lett. 99, 4, 1–6, doi: 10.1209/02955075/99/44004. [CrossRef] [Google Scholar]
 Sinha S. Bender A.T., Danczyk M., Keepseagle K., Prather C.A., Bray J.M., Thrane L.W., Seymour J.D., Codd S.I., Hansen A. (2017) Effective rheology of twophase flow in threedimensional porous media: experiment and simulation, Transp. Porous Media, doi: 10.1007/s1124201708744. [PubMed] [Google Scholar]
 Taber J.J., Martin F.D., Seright R.S. (1997) EOR screening criteria revisited − part1: Introduction to screening criteria and enhanced recovery field projects, SPE Reserv. Eng. 189–198 SPE35385. [CrossRef] [Google Scholar]
 Taber J.J., Martin F.D., Seright R.S. (1997) EOR screening criteria revisited − part2: Applications and impact of oil prices, SPE Reserv. Eng. 199–205 SPE39234. [CrossRef] [Google Scholar]
 Tallakstad K.T., Knudsen H.A., Ramstad T., Løvoll G., Maløy K.J., Toussaint R., Flekkøy E.G. (2009) Steadystate twophase flow in porous media: statistics and transport properties, Phys. Review Lett. 102, 074502, 1–4, doi: 10.1103/PhysRevLett.102.074502. [CrossRef] [PubMed] [Google Scholar]
 Tsakiroglou C.D., Avraam D.G., Payatakes A.C. (2007) Transient and steadystate relative permeabilities from twophase flow experiments in planar pore networks, Adv. Water Resour. 30, 1981–1992, doi: 10.1016/j.advwatres.2007.04.002. [CrossRef] [Google Scholar]
 Tsakiroglou C.D., Aggelopoulos C.A., Terzi K., Avraam D.G., Valavanides M.S. (2015) Steadystate twophase relative permeability functions of porous media: A revisit, Int. J. Multiph. Flow 73, 34–42, doi: 10.1016/j.ijmultiphaseflow.2015.03.001. [CrossRef] [Google Scholar]
 Tzimas G.C., Matsuura T., Avraam D.G., van der Brugghen W., Constantinides G.N., Payatakes, A.C. (1997) The combined effect of the viscosity ratio and the wettability during forced imbibition through nonplanar porous media, J. Colloid Interface Sci. 189, 27–36, doi: 10.1006/jcis.1996.4658. [CrossRef] [Google Scholar]
 Valavanides M.S. (1998) Macroscopic theory of twophase flow in porous media based on integration of pore scale phenomena, PhD Dissertation, University of Patras, National Archive of PhD Theses − National Documentation Center, doi: 10.12681/eadd/11044. [Google Scholar]
 Valavanides M.S., Constantinides G.N., Payatakes A.C. (1998) Mechanistic model of steadystate twophase flow in porous media based on ganglion dynamics, Transp. Porous Media 30, 267–299, doi: 10.1023/A:1006558121674. [CrossRef] [Google Scholar]
 Valavanides M.S., Payatakes A.C. (2000) A truetomechanism model of steadystate twophase flow in porous media, including the contribution of the motion of ganglia and droplets, in: Bentley L.R. et al. (eds.), Computational Methods in Water Resources XIII, 1, ISBN 9058091236, A.A. Balkema, Rotterdam, The Netherlands, pp. 239–243. [Google Scholar]
 Valavanides M.S., Payatakes A.C. (2001) Truetomechanism model of steadystate twophase flow in porous media using decomposition into prototype flows, Adv. Water Resour. 24, 3–4, 385–407. [CrossRef] [Google Scholar]
 Valavanides M.S., Payatakes A.C. (2002) Effects of Pore Network Characteristics on SteadyState TwoPhase Flow Based on a TruetoMechanism Model (DeProF), SPE78516MS, 10th ADIPEC Abu Dhabi International Petroleum Exhibition & Conference, Abu Dhabi, United Arab Emirates, October 13–16, pp. 379–387, https://doi.org/10.2118/78516MS. [Google Scholar]
 Valavanides M.S., Payatakes A.C. (2003) Prediction of optimum operating conditions for steadystate twophase flow in pore network systems using the DeProF truetomechanism theoretical model, SCA200318, 2003 Int. Symp. Soc. Core Anal., Pau, France, 21–25 September. [Google Scholar]
 Valavanides M.S., Payatakes A.C. (2004) Wetting Film Effects on SteadyState TwoPhase Flow in Pore Networks using the DeProF Theoretical Model, SPE88713MS, 11th ADIPEC Abu Dhabi International Petroleum Exhibition & Conference, Abu Dhabi, United Arab Emirates, October 10–13, https://doi.org/10.2118/88713MS. [Google Scholar]
 Valavanides M.S. (2012) Steadystate twophase flow in porous media: review of progress in the development of the DeProF theory bridging pore − to statistical thermodynamics − scales, Oil Gas Sci. Technol. 67, 5, 787–804, doi: 10.2516/ogst/2012056. [CrossRef] [Google Scholar]
 Valavanides M.S. (2014) Operational efficiency map and flow characterization for steadystate twophase flows in porous media, Soc. Core Anal. Symp. − SCA2014, Avignon, France, Sept. 814, http://users.teiath.gr/marval/publ/Valavanides_SCA2014047.pdf. [Google Scholar]
 Valavanides M.S. (2018) Review of steady state twophase flow in porous media: independent variables, universal energy efficiency map, critical flow conditions, effective characterization of flow and pore network. Transp. in Porous Media online, 155, doi: 10.1007/s1124201810261]. [Google Scholar]
 Valavanides M.S., Daras T. (2016) Definition and counting of configurational microstates in steadystate twophase flows in pore networks, Entropy 18, 054, 1–28, doi: 10.3390/e18020054. [CrossRef] [Google Scholar]
 Valavanides M.S., Totaj E., Tsokopoulos M. (2016) Energy efficiency characteristics in steadystate relative permeability diagrams of twophase flows in porous media, J. Petrol. Sci. Eng. 147, 181–201, doi: 10.1016/j.petrol.2016.04.039. [CrossRef] [Google Scholar]
 Van de Merwe W., Nicol W. (2009) Trickle flow hydrodynamic multiplicity: Experimental observations and porescale capillary mechanism, Chem. Eng. Sci. 64, 1267–1284, doi:10.1016/j.ces.2008.10.069. [CrossRef] [Google Scholar]
 Vizika O., Avraam D.G., Payatakes A.C. (1994) On the role of viscosity ratio during lowcapillarynumber forced imbibition in porous media, J. Colloid Interface Sci. 165, 386–401. [CrossRef] [Google Scholar]
 Youssef S., Rosenberg E., Deschamps H., Oughanem R., Maire E., Mokso R. (2014) Oil ganglia dynamics in natural porous media during surfactant flooding captured by ultrafast xray microtomography, SCA201423, 2014, Int. Symp. Soc. Core Anal., Avignon, France, 11–18 September. [Google Scholar]
All Tables
Occurrence probabilities, f_{j}, and respective reduced size diameters for the 5 classes of chambers, D_{C}_{j}, and throats, D_{T}_{j}, of the 3D pore network used in the DeProF simulations. All dimensions are normalized to the lattice constant, ℓ = 1221 μm.
Volume fraction of the connected pathway flow (CPF) domain, β, against fractional flow of oil through it, η_{o,C}, for various flow conditions, (Ca, r), and viscosity ratios κ.
Reduced pressure gradient values, x, for different flow conditions, (Ca, r), and viscosity ratios, κ.
All Figures
Figure 1 (a) Schematic representation of the actual flow (left sketch) and its theoretical decomposition into prototype flows: Connectedoil Pathway Flow (CPF) and DisconnectedOil Flow (DOF) (right sketch), (b) A microscopic scale representation (snapshot) of a DOF region. An oil ganglion of size class 5 is shown. For simplicity, all cells are shown identical and the lattice constant is shown expanded (chambers and throats have prescribed size distributions). The dashed, light grid lines define the network unitcells. The thick dashed line separates the Ganglion Dynamics (GD) cells domain and the DropTraffic Flow (DTF) cells domain. 

In the text 
Figure 2 The model of the ganglion size distribution expressed through equation (15) by three parameters, 

In the text 
Figure 3 Typical domains of physically admissible flow configurations, as predicted by DeProF model simulations. For any fixed, externally imposed, flow conditions, (Ca, r), the twophase flow visits a canonical ensemble of physically admissible flow configurations depicted by the cloud of small, red balls. Here, Ca = 1.19 × 10^{−6}, κ = 1.45 and flow rate ratio values, r, span 2 orders of magnitude. A unique triplet of values (S_{w}, β, ω), the large, black ball, is obtained by averaging over the ensemble and defines the mass center of the red cloud. Note that the shape and position of the red cloud change with flow conditions. 

In the text 
Figure 4 Input/output block diagram in the DeProF model implementation. The variables describing the structure of the interstitial flow are presented in Section 3. 

In the text 
Figure 5 Typical jikclass unit cell of the model pore network used in the simulations. The thick inclined arrow represents the orientation of the lattice skeleton relative to the macroscopic pressure gradient. 

In the text 
Figure 6 (a) A typical 3D cubic pore network. The lattice constant, ℓ, compared to the size of the throats and chambers, is shown expanded. The thick arrow indicates the direction of the macroscopic flow, (b) The pore network cubic lattice with one main diagonal parallel to the macroscopic flow direction. The two parallel equilateral triangles are perpendicular to the macroscopic flow, (c) Equilateral triangles align to form a checkerboard pattern and exactly fill any frontal area perpendicular to the macroscopic flow. 

In the text 
Figure 7 A closeup sketch view of Figure 1(b) indicating the different types of unit cells (u.c.) in the DOF domain. The classification of the ganglion unit cells (‘E,G’, ‘X,G’ and ‘C,G’) is according to their conductance state. 

In the text 
Figure 8 (a) A ganglion of size class 4 at a least expected to get mobilized configuration in the time interval between two rheons. (Emphasis is on the endgate unit cells), (b) The two configurations for which an oil droplet may become stranded. In both cases, the driving force for mobilization is from left to right. represents the idealized radii of curvature of the o/w interfaces for advancing and receding contact angles. 

In the text 
Figure 9 Mean mobilization probabilities for various types and sizes of oil blobs, as members of a dense population in terms of the reduced macroscopic pressure gradient, x. Continuous red line shows the mobilization probability of droplets, m^{D}, discontinuous black lines the mobilization probability of ganglia of various size classes, . 

In the text 
Figure 10 Expected, ensemble average values of the reduced pressure gradient, x, for different values of the capillary number, Ca, flow rate ratio, r, and viscosity ratio , Arrows indicate increasing Ca values: →10^{−8} ≤ Ca ≤ 10^{−4} and    > 4 × 10^{−7} ≤ Ca ≤ 10^{−4}. 

In the text 
Figure 12 Expected, ensemble average values of the flow arrangement variables for different values of the capillary number, Ca, and viscosity ratio, κ. (a) Saturation, S_{w}; (b) volume fraction of connected pathway flow (CPF), β; (c) volume fraction of ganglion dynamics flow (GD) within the disconnected oil flow (DOF), ω. Arrows indicate increasing Ca values: →10^{−8} ≤ Ca ≤ 10^{−4} and    > 4 × 10^{−7} ≤ Ca ≤ 10^{−4}. 

In the text 
Figure 13 Expected, ensemble average values of the reduced mechanical power dissipation, W, and its partition into the individual interstitial, constituent flows, for different values of the capillary number, Ca, flow rate ratio, r, and viscosity ratio, . Arrows indicate increasing Ca values: →10^{−8} ≤ Ca ≤ 10^{−4} and    > 4 × 10^{−7} ≤ Ca ≤ 10^{−4}. 

In the text 
Figure 14 Expected, ensemble average values of the mean ganglion sizes for different values of the capillary number, Ca, and viscosity ratio, κ. a) Mean size of all ganglia (mobilized + stranded), ; b) mean size of mobilized ganglia, . Arrows indicate increasing Ca values: →10^{−8} ≤ Ca ≤ 10^{−4} and    > 4 × 10^{−7} ≤ Ca ≤ 10^{−4}. 

In the text 
Figure 15 Expected, ensemble average values of the reduced o/w interface surface area for different values of the capillary number, Ca, flow rate ratio, r, and viscosity ratio, . (1st row) ganglion dynamics, A_{ow,G}; (2nd row) oil droplets, A_{ow,D}; (3rd row) Disconnected oil flow, A_{ow}; (4th row) overall. Arrows indicate increasing Ca values: →10^{−8} ≤ Ca ≤ 10^{−4} and    > 4 × 10^{−7} ≤ Ca ≤ 10^{−4}. 

In the text 
Figure 16 Expected, ensemble average values of the coefficient (fraction) of oil fragmentation, f_{OF}, (or interfacial saturation) for different values of the capillary number, Ca, flow rate ratio, r, and viscosity ratio, . Arrows indicate increasing Ca values: →10^{−8} ≤ Ca ≤ 10^{−4} and    > 4 × 10^{−7} ≤ Ca ≤ 10^{−4}. 

In the text 
Figure 17 Expected, ensemble average macroscopic fractional flows through connected pathways, η_{o,C}, ganglion dynamics, η_{o,G}, and drop traffic, η_{o,D}, flows. Arrows indicate increasing Ca values: →10^{−8} ≤ Ca ≤ 10^{−4} and    > 4 × 10^{−7} ≤ Ca ≤ 10^{−4}. 

In the text 
Figure 18 Expected, ensemble average values of the o/w surface transport variables for different values of the capillary number, Ca, flow rate ratio, r, and viscosity ratio, . a) Reduced superficial velocity of o/w surface, U_{ow,DOF}, b) Fraction of o/w surface transport through ganglion dynamics , c) Fraction of o/w surface transport through drop traffic flow,, Arrows indicate increasing Ca values: →10^{−8} ≤ Ca ≤ 10^{−4} and    > 4 × 10^{−7} ≤ Ca ≤ 10^{−4}. 

In the text 
Figure 19 Effect of the flow rate ratio on the energy efficiency of the process (oil output per kW of power dissipated by the process), f_{EU} = r/W, for different values of the capillary number, Ca, flow rate ratio, r, and viscosity ratio, . Arrows indicate increasing Ca values: →10^{−8} ≤ Ca ≤ 10^{−4} and    > 4 × 10^{−7} ≤ Ca ≤ 10^{−4}. 

In the text 
Figure 11 Diagrams of the ganglion size distributions, pertaining to moderate physically admissible flows for different macroscopic flow conditions, (a1), (b1), (c1), and to extreme physically admissible flows for the same macroscopic flow conditions (a2), (b2), (c2). Tabulated data indicate the values of flow arrangement variables, and of corresponding ganglion size distributions parameters, . The values of flow conditions, (Ca, r), and average flow arrangement variables,(S_{w}, β, ω), are provided in the upper two rows of each tablet. 

In the text 
Figure 20 The set of critical flow conditions, , detected by the DeProF simulations, for which the energy efficiency takes maximum values (ref. Fig. 19). 

In the text 