Oil Fragmentation, Interfacial Surface Transport and Flow Structure Maps for Two-Phase Flow in Model Pore Networks. Predictions Based on Extensive, DeProF Model Simulations

— In general, macroscopic two-phase ﬂ ows in porous media form mixtures of connected-and disconnected-oil ﬂ ows. The latter are classi ﬁ ed as oil ganglion dynamics and drop traf ﬁ c ﬂ ow, depending on the characteristic size of the constituent ﬂ uidic elements of the non-wetting phase, namely, ganglia and droplets. These ﬂ ow modes have been systematically observed during ﬂ ow within model pore networks as well as real porous media. Depending on the ﬂ ow conditions and on the physicochemical, size and network con ﬁ guration of the system ( ﬂ uids and porous medium), these ﬂ ow modes occupy different volume fractions of the pore network. Extensive simulations implementing the DeProF mechanistic model for steady-state, one-dimensional, immiscible two-phase ﬂ ow in typical 3D model pore networks have been carried out to derive maps describing the dependence of the ﬂ ow structure on capillary number, Ca , and ﬂ ow rate ratio, r . The model is based on the concept of decomposition into prototype ﬂ ows. Implementation of the DeProF algorithm, predicts key bulk and interfacial physical quantities, fully describing the interstitial ﬂ ow structure: ganglion size and ganglion velocity distributions, fractions of mobilized/stranded oil, speci ﬁ c 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 non-wetting phase in disconnected form is signi


NOTATION ∼: A tilde embellishment indicates a dimensional variable
No tilde embellishment denotes a dimensionless variable F 0 : A prime indicates a physically admissible value of variable F (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 Conductance of the various unit cells in the GD and DTF domains P 1.3

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 geo-sequestration, soil drying/wetting or pollution remediation, operation of packed-bed reactors, operation of proton exchange membrane fuel cells, etc. Immiscible two-phase 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 non-wetting phase (oil) produces small fluidic elements, called ganglia and droplets.Ganglia are micro-fluidic 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 non-wetting and the solid phase Porous Medium (p.m.) create interfaces separating the non-wetting fluidic elements from the wetting phase.The collective effect of inherent flow mechanisms, prevailing over the pore (micro-to-millimeter) scale domain and described by Stokes and Young-Laplace equations, results in the peculiar phenomenology of two-phase flows in porous media observed on the core, fracture, or field/production scales.
In general, macroscopic two-phase flows in porous media form mixtures of connected-and disconnected-oil 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 andPayatakes, 1996 and2000;Knudsen et al., 2002;Knudsen and Hansen, 2002;Al-Gharbi 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 disconnected-oil 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 oil-flow, 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 oil-water-p.m. system (Avraam et al., 1994;Avraam andPayatakes, 1995 and1999;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 non-linear, and to be strongly affected by the physical parameters that pertain to the fluid-fluid interfaces.
Focusing on steady-state, immiscible, two-phase 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 Connected-Oil Pathway Flow (CPF) and Disconnected-Oil 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 steady-state twophase flows in pore networks, Valavanides andPayatakes (2000 and2001;also Valavanides, 1998 and2012) developed the DeProF mechanistic model.The model is based on the concept of decomposition in prototype flows (see next section).It accounts the pore-scale mechanisms and the network-wide cooperative effects.The sources of non-linearity, 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 pore-to-macro 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 two-phase flow in pore networks, to understand the role of disconnected flow and to reveal if any dynamic similarities apply across different oil-water-p.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 carried-out 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 two-phase flows in a typical 3D model pore network of the chamber-and-throat 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 two-phase flows in porous media in mixtures of connected and disconnected oil flow (connected-oil 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 non-wetting 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 non-wetting 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 non-wetting 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 nano-particles (Ryoo et al., 2010) and other species (Serres-Piole et al., 2012), or the rate of dissolution of the non-wetting 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 trickle-bed 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.

THE DEPROF HIERARCHICAL MECHANISTIC MODEL FOR STEADY-STATE TWO-PHASE 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 steady-state two-phase flow in pore networks, by predicting the reduced macroscopic pressure gradient, x, given the values of physical quantities defining the flow conditions and the oil-water-p.m. system properties.

The independent variables for steady-state twophase flow in pore networks
The conditions of steady-state two-phase flow in pore networks are defined by the values of the flow rate intensities, i.e. superficial velocities, Ũ o and Ũ w , conventionally defined as oil and water flow rates, qo and qw , divided by the p.m. cross-sectional area, Ã.The set of superficial velocities may be appropriately reduced and replaced by a set of dimensionless variables, comprising the capillary number, and the oil/water flow rate ratio, where mw is the viscosity of water and g ow is the oil/water interfacial tension.
The oil-water-p.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., k; the oil/water viscosity ratio; the dynamic advancing and receding contact angles of the o/w interface on the pore walls, u A and u 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, chamber-to-throat size correlation factors, etc).

The concept of Decomposition into Prototype Flows (DeProF)
In the general case of steady-state, one-dimensional, 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, z, with constant flow rates, qo and qw , 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.
Water is considered to be the wetting phase and always retains its connectivity, while oil (as the non-wetting phase) gets disconnected to form fluidic elements of various sizes: from infinite-long and slim connected-oil 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, -Connected-oil Pathway Flow (CPF) and; -Disconnected Oil Flow (DOF).
The latter can be further decomposed in -Ganglion Dynamics (GD) flow and; -Drop Traffic Flow (DTF).
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 one-phase flow.The p.m. volume fraction occupied by the connected oil is denoted by b.The Disconnected-Oil Flow (DOF) regime is defined as the region composed of the rest of the pore network, so the associated volume fraction equals (1 À b).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 re-accounted 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 non-wetting phase, the REV pertaining to immiscible two-phase 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 v, and is called the GD domain fraction.The DTF domain fraction in the DOF region equals (1 À v).
S w , b and v 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 , b and v that conform with the externally imposed flow conditions (Ca, r).
The flow analysis is carried-out 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 micro-to-macroscopic scale consistency relations-see next subsection.
Scale-wise, 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 (scaled-up) values of pore-scale (P) variables, (e.g. the statistical average velocity of an i-class 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.

The DeProF model equations
Using the DeProF model, one can obtain the solution to the problem of steady-state two-phase flow in porous media in the form of the following transfer function where, is the reduced macroscopic pressure gradient, i.e. the actual pressure gradient divided by the pressure gradient of an equivalent one-phase flow of water at superficial velocity equal to Ũ w [the 2nd component of the product in Equation ( 5)].
Referring to Figure 1, the externally imposed flow rates of oil, qo , and water, qw , are partitioned in the prototype flows.In particular, qo;C , qo;DOF are the mean oil flow rates in CPF and DOF and qw;DOF the mean water flow rate in DOF (as CPF is oil-saturated, qw;DOF ¼ qw ).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 U o;C ¼ qo;C =q o , and the reduced oil and water mean superficial velocities defined by U m;DOF ¼ qm;DOF =q m , where m = o, w.
The flow analysis produces the following system of equations (Valavanides and Payatakes, 2000;2001) Uo 1 2 where n G i ; i ¼ 1; :::; I max is the reduced number density of the i-class ganglia and I max is the upper truncation ganglion size class.N G i is the number of unit cells occupied by any iclass ganglion.u G i ð¼ ũG i = Ũ w Þ is the reduced mean ganglion velocity i-class ganglia.m G i is the mean mobilization probability of the i-class ganglia.P b a;i and P D a are the occurrence probabilities of the various unit cells in the GD and DTF domains respectively.g b a;i and g D a are the reduced conductance of the various unit cells in the GD and DTF domains respectively.
Note that m D , u G i , m G i , g b a;i and g D a depend on x, as described by semi-analytical expressions of momentum balance derived over the pore scale.P b a;i , P D a , g b a;i and g D a also depend on the reduced pressure gradient, x, as well as on the unit cell geometry (index 'a'), the two-phase flow pattern in it (index 'b') and the ganglion size (index 'i').P D a and g D a 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 Darcy-type 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 ganglia-occupied 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 freedomby 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 one-phase flow in the DOF (GD and DTF) region and implicitly represents the transfer function for this region.
The above system of 8 equations relates ( , along with the three FAV, namely S w , b, v.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 semi-analytically and/or numerically infused in (13).
Note that the most difficult part in implementing the DeProF model is in estimating the dependence of P b a;i , P D a , g b a;i and g D a on the given pore network geometry.This may be accomplished either by furnishing appropriate analytical expressions (if possible) for P b a;i , P D a , g b a;i and g D a or by delivering appropriate maps implementing advanced numerical techniques (CFD or L-B) and/or sophisticated laboratory experiments.Analytical expressions have already been derived in the past for the structure/geometry of the chamber-and-throat type model pore network described in Section 2 (Valavanides, 1998;Valavanides and Payatakes, 2001).These have been used in the present work.

The Discrete Fluidic Elements of the Non Wetting Phase
Discrete fluidic elements are formed by the disconnection and immiscible dispersion of the non-wetting 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 non-wetting 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 infinite-many pore unit-cells.
In general, three classes of discrete fluidic elements are considered, depending on their size: connected-oil pathways comprising the Connected Pathway prototype Flow (CPF), oil ganglia, occupying one-to-many unit-cells comprising the GD prototype flow, and oil droplets dispersed within water saturated unit-cells 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 Connected-Oil Pathways (COP), i.e. connected-oil unit-cells with infinite length.COPs may occupy or arrange themselves within any b region of the pore network and remain aligned to the macroscopic flow.Of a total number of M unit-cells in a reference volume, bM unit-cells are occupied by connected oil (the ensemble of Connected-Oil Pathways) and allow the rest (1 À b) M of the unit-cells to host the Disconnected-Oil Flow (DOF) unit-cells, Figure 1a.
DOF occupies the remaining unit-cells within any (1 À b) 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 oil-droplets and water.Oil ganglia unit-cells take up a volume fraction of (1 À b) v of the pore network, leaving to the DTF cells a complementary volume fraction of (1 À b) (1 À v).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 unit-cells [Figs.1b and 2] are not indistinguishable in the sense that each ganglion unit-cell is part of a certain ganglion size class.Oil ganglia have a sharply decreasing volume size distribution, defined as where N d is the network dimensionality (=3 for 3D networks), Ṽ uc is the average volume of unit cells and n G i denotes the ratio of the total number of i-class 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 where n 1 , n 2 , I max and z are parametric values that define the particular physically admissible flow configuration.I max is the maximum attainable size of a ganglion and 0 < z < 1 provides a measure of the sharpness of the ganglion size distribution, see Figure 2. 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 break-up 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., 2014;Youssef et al., 2014;Georgiadis et al., 2013;Datta et al., 2014;Oughanem et al., 2015;Armstrong et al., 2016).The model of the ganglion size distribution expressed through Equation ( 15) by three parameters, n G 1 ; n G 2 ; z: Indicative DeProF model predictions of physically admissible ganglion size distributions, ðn G 2 ; n G 2 ; zÞ, are presented in the diagrams in Section 3.1, along with the corresponding values (tabulated) of the FAV (S w , b, v).The expected, average values of the flow arrangement variables, (S w , b, v), for different values of the capillary number, the flow rate ratio and the viscosity ratio, are presented in the diagrams in Section 3.2.

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 , b, v} 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 , b, v), 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 , b, v} 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 sub-domain of the {S w , b, v} space, denoted V 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.The measure ("volume") of V PAS is given by where V PAS (Ca, r) stands for integration carried over the physically admissible ranges in {S w , b, v} for the imposed pair of (Ca, r) values.The volume of V 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, V PAS , a unique solution for the macroscopic flow is obtained.For any quantity, F ', the corresponding expected mean macroscopic flow quantity, F is defined as A prime is used to denote physically admissible values of any quantity.The symbol without a prime is reserved for the expected, ensemble-averaged value of that quantity.
As with any macroscopic physical quantity, unique triplet of values for (S w , b, v) can be obtained by averaging over the domain of physically admissible solutions, using Equation ( 17), with F replaced by S w , b and v.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.
In computational practice, the size of the domain of physically admissible solutions, V V PAS , as well as any expected mean macroscopic flow quantity, F, 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 Â b Â v} pertaining to the FAV, is partitioned in corresponding N S Â N b Â N v 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 ðS 0 w ; b 0 ; v 0 Þ i and ðS 0 w ; b 0 ; v 0 Þ iþ1 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 picked-up 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 re-discretizing 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).
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 cause-effect 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 non-wetting phase will settle-in to comply with the externally imposed flow conditions, ð Ũ o ; Ũ w Þ 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).Additional comments are provided in Section 4.9 'Discussion of results' in the present work.

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.

Pore Network Geometry
The simulations are performed on the three-dimensional (3D) cubic pore network models of the chamber-and-throat 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 andPayatakes (1995, 1999).
The 3D cubic lattice pore network is of the "ball-andstick" 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 jik-class unit cell comprises two chambers of classes j and k with diameters D Cj and D Ck , interconnected with throat-i, 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 dashed-dotted lines are the two of the three octahedron diagonals (Fig. 5).
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 s 3D = 6.For the particular network dimensions, the volume porosity is e = 0.0464 and its absolute permeability is k ¼ 8:965mm 2 .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 unit-cell occupies an equilateral octahedron of pore network space with diagonal lengths equal to ℓ, edge lengths equal to ℓ ffiffi ffi 2 p =2 and volume ℓ 3 =6.The geometry of the pore space within every unit-cell is defined by the shape of the chambers and throats (Fig. 5).For this particular type of cubic lattice networks, each unit-cell 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 cross-sectional 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 fill-up any cross-section perpendicular to the macroscopic flow (frontal area).Focusing on the macroscopic flow direction in Figure 6, in every unit-lattice there exist two parallel equilateral triangles, oriented perpendicular to the main diagonal and to the macroscopic flow.These are depicted in Figure 6 and its frontal area is equal to ffiffi ffi 3 p ℓ 2 =4.An important parameter associated with the 'zigzag' motion of the ganglion is the tortuosity of its spine, x 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,

Types and number of unit cells in the Disconnected Oil Flow domain (DOF)
Ganglia belonging in the same size-class 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 drop-traffic flow cells are identical (momentum-wise) since they contain water and uniformly dispersed oil droplets.

Ganglion unit cells
The number N G i of unit cells occupied by each i-class ganglia are given by 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' (end- gate, ganglion), see Figure 7.The only potential end-gate unit cells are those located at the foremost part (for the hygron) and at the aft-most 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 end-gate cells nor core cells) contain oil and water separated by virtually immobile menisci and are called side-gate unit cells.Ganglion migration induces a corresponding flow rate, qG uc;i , in the core cells and an equal flow rate in the end-gate cells.This flow must be consistent with the mean velocity of the mass center of the i-class ganglion, ũG i , i.e., where x G is the network tortuosity, defined in the previous subsection.
The rest of the cells of a mobilized ganglion (side-gate unit cells) are considered as non-conducting, and they are denoted with the index 'X,G' ('non-conducting, ganglion').All of the N G i cells of any i-class stranded oil ganglion are 'X,G' cells.
The fractional distributions of the various flow configurations of any i-class 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, Non-conducting ganglion cells (note here, that this fraction takes into account all the cells of stranded ganglia), End ganglion cells À conducting, through which hygrons take place, Now, as expected, 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 micro-roughness.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 one-to-one correspondence with pore/throat size classes.Nevertheless, as the population of oil droplets migrates downstream, they are mixed and there is no correlation between pore-size and resident dropletsizes.The flow set-up in the DTF domain unit-cells is considered similar to oil-in-water 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.

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 log-normal 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.
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 summing-up 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 low-end 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.

VARIABLES DESCRIBING THE INTERSTITIAL FLOW STRUCTURE
An inherent characteristic of immiscible two-phase flow in porous media, is the disconnection of the non-wetting 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, break-up and pinch-off of the non-wetting phase.These phenomena are counterbalanced À to a certain extent À by collision and eventual coalescence of non-wetting blobs (oil ganglia and oil droplets).The combined effect of these events leads to the collective behavior of the disconnected-oil 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 À scale-up.
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 disconnected-oil 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 bottom-up 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 summed-up over the various flow configurations and/or the blob populations (ganglia and droplets) or the entire connected or disconnected flow domains.

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, n 0 G i ; i ¼ 1; :::; I 0 max , 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, V 0 G m , is estimated as the ratio of the entire oil ganglia mass over the entire ganglia population, as where V G i is the volume of oil comprising any ganglion of class-i or, equivalently, the volume of i-class ganglia unit cells.
Similarly, the mean ganglion size pertaining to the mobilized ganglia population, V 0 MG m , is estimated as the ratio of the entire oil mass of mobilized ganglia over the population of mobilized ganglia, as Consequently, the mean ganglion size pertaining to the population of stranded ganglia, V 0 SG m , is estimated as For simplicity, such estimates are not presented here.

Specific (p.u.v.p.m.) surface areas of o/w-interface 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/w-interface gets displaced across the front-most gate unit cell hosting the head (bow) of the ganglion (Fig. 7).During a hygron, an advancing o/w interface gets displaced across the aft-most 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 i-class ganglion is denoted as ÃmMG ow;i .This mobile o/w interface comprises exclusively the respective receding and advancing o/w interfaces in the two gate unit cells (E,G).The mean-across-the gate-u.c.surface of o/w-interfaces, A rec ow and A adv ow , are computed for all jik-classes of unit cells (combinations of j and k chambers connected with i-throats) during the initial preprocessing steps (in the DeProF simulations) and considering the respective À advancing and receding À contact angles.
The rest of the o/w-interfaces 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, Mean surface area of the stranded o/w interfaces in the X, G u.c. of an i-class ganglion, whereby the mean is considered over all X,G unit cell classes, Total mean surface are of the o/w interface of an i-class mobilized ganglion Total mean surface area of the o/w interface of an i-class stranded ganglion, whereby the mean is considered over all unit cell classes, In the previous expressions, Ãchamber ow ; Ãthroat ow , 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), A In all relevant expressions, the index notation consistently determines the physical quantities as follows: Gganglia, D -Droplets, (ow)MMG/SMGmobilized/stranded interface of mobilized ganglia, (ow)MD/SDmobilized/stranded interface of droplets, (ow)Gtotal interface area of ganglia etc.In particular, the mobilized o/w surface area of mobilized ganglia p.u.v. of GCD is estimated by the stranded o/w surface area of mobilized ganglia p.u.v. of GCD is the total o/w surface area of mobilized ganglia p.u.v. of GCD is Similarly, the o/w surface area of stranded ganglia p.u.v. of GCD, is The stranded o/w surface area of ganglia (mobilized and stranded) p.u.v. of GCD, is The total o/w surface area of ganglia (mobilized and stranded) p.u.v. of GCD, is Specific (p.u.v.p.m.) surface areas of o/w-interface in the DTF region The average (mean) oil drop volume is The average number of oil drops p.u.v. of DTF domain is The o/w surface area of mobilized drops p.u.v. of DTF domain is while the o/w surface area of stranded drops p.u.v. of DTF domain is where A mD ow and A sD ow 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, A mD ow , is computed as drops attain all possible configurations during their migration.Every n-class drop is screened and then its mean o/w surface area is computed as it crosses a jik-class unit cell.Similarly, the mean o/w surface area of stranded drops, A sD ow , is computed as drops attain all possible configurations during their stranding.Every n-class drop is screened to compute its mean surface of the o/w surface area as it blocks a jik-class 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 The total o/w surface area p.u.v. of DOF domain is The total o/w surface area p.u.v. of the whole pore network is 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, f 0 OF , 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 one-to-one correspondence with throat size classes.The minimum size of an oil-droplet is, whereas the corresponding o/w surface area is 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, Ṽ D;1 ℓ 3 , The respective specific o/w surface area, Ã0 ow;crit , would comprise spherical surfaces pertaining to the mimimum size drops, of diameter equal to the diameter of the minimum throat, Dh1 , therefore, The coefficient (or, fraction) of oil fragmentation is defined as It is straightforward to derive its limiting values as 0 f 0 OF 1, and based on its definition, the coefficient of oil fragmentation, f OF , could also be called 'interfacial saturation'.Expected, ensemble average values of f 0 OF are presented in Section 4.

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 disconnected-oil 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, where U 0 o;C 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, h o,C , GD, h o,G , and DTF, h 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, Ũ 0 ow;G , may be estimated as the total momentum of ganglia o/w surface divided by the total o/w surface area, and then, reduced by the macroscopic velocity of oil, it becomes: 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 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, U 0 o;D .

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 Fractional o/w surface transport through GD and DTF The fractional transport of the o/w interfaces through GD and DTF, j 0 ow;G and j 0 ow;D , may be estimated as the momentum of ganglia o/w surface divided by the total o/w surface area,

SIMULATIONS OF STEADY-STATE 2-PH FLOW IN PORE NETWORKS
A large number of simulations of steady-state two-phase flows in the three-dimensional (3D) chamber-and-throat 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, mw ¼ 0:94 Â 10 À3 Pas, oil/ water interfacial tension g ow ¼ 25 Â 10 À3 N=m, advancing (A) and receding (R) oil-water/p.m. contact angles, u 0 A ¼ 43deg and u 0 R ¼ 37deg.Five typical systems were examined: favorable, k = 0.33, 0.67, indifferent k = 1.00 and unfavorable, k = 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. steady-state two-phase flow is not possible/sustainable.Because of space limitations, albeit 5 different viscosity ratio values were examined, only the results pertaining to k = 0.33, k = 1.00 and k = 3.00 are presented.The flow structure maps of the other two systems (k = 0.67 and 1.50) fall in between.
The aforementioned diagrams (Figs. 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 iso-capillary curves shift from right to left and extend over narrower logr domains.The iso-capillary 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 iso-capillary 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 (k = 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.8provide 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.

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, k.
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 two-phase 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 steady-state does not scale linearly with the flow rate.
The DeProF model predictions in Figure 10 are consistent with recent experimental studies on steady-state two-phase flows in glass beads in Hele-Shaw 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 glass-etched pore network models and in sand-pack column experiments (Tsakiroglou et al., 2015), observing that the non-linear 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 are 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.

Reduced mechanical power dissipation
The reduced rate of mechanical energy dissipation within the process, W, is defined as (Valavanides et al., 2016) where is the mechanical power dissipation for the equivalent one-phase flow of water with superficial velocity, Ũ w .The mechanical power dissipated in Connected Pathway Flow (CPF), Ganglion Dynamics (GD) and Drop-Traffic Flow (DTF) are given respectively by 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, k , are presented in the diagrams of Figure 13.

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, k, are presented in the diagrams of Figure 14.The mean size pertaining to the entire ganglia population (mobilized and stranded), V G m , is presented in (a), whereas the mean size of the mobilized ganglia,V mG m , in (b).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, k, 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 connected-oil pathways (no o/w interface).This drop becomes more intense by shifting into higher iso-capillary curves, as observed by the steeper drop from higher values in f OF as flow conditions shift to higher iso-capillary 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.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 Disconnected-Oil Flow (DOF), the contribution of DTF is generally more intense than that through GD.At relatively low-Ca/low-r 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 takes-over a substantial fraction of the oil flow over droplets and ganglia.Then, with increasing Ca, this take-over 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, j ow,D , and GD, j 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.

Energy efficiency of steady-state two-phase flow in porous media processes
The efficiency of any steady-state two-phase 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] as 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 oilwater-p.m. system (Valavanides, 2014 andValavanides, 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, r Ã i .The sets Ca i ; r Ã i È É 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, dash-dotted line in red indicates the theoretically predicted values of the (global) maximum attainable energy efficiency (the o/w/pm system's 'efficiency ceiling'), .These conditions are theoretically met for 'infinite' Ca, in which case, o/w interfaces have no effect at all to the viscosity dominated flow.
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 steady-state two-phase 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 trading-off between CPF, GD and DTF and self-adjusting the connected versus disconnected moving-oil 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 steady-state two-phase flows in porous media-179 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 two-phase flows in real 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.

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 steady-state two-phase 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 saturated-flow 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), b !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. cross-section, as is indicated by the rather low v and high (1 À v) 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 non-existing [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. cross-section, 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 V mG m ≈ 10 [diagrams in row (b)], large enough, if compared to the mean size of all ganglia,V G m ≈ 1:5 [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 two-phase flow (À2 < logr < 0), from pure DOF, (1 À b) ! 1, and borderline CPF b !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. cross-section À over GD, as is indicated by the low v and high (1 À v) 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 k = 0.33), its contribution in the oil transport is dramatic as it suddenly steps-up 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 k = 0.33.With further increase in r the fractional flow of oil through CPF climbs-up towards 100%, until two-phase 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 V G m ≈ 1:0 À Á , see Figure 14a.The mean size of mobilized ganglia also small 2 < V mG m < 3 À Á , 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, h o,C , are larger than the corresponding values of the volume fractions, b.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 connected-oil pathways, b, Figure 12b, with those of the fractional flow of oil through CPF, h o,C , Figure 17b.Indicative values are tabulated for typical flow conditions and viscosity ratios, see Table 2.The 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 two-fluids have equal viscosities, two-phase 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 low-Ca 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 high-Ca regime, mobilization is more probable, therefore two-phase flow is possible for low-r values.In that case (k = 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 non-wetting 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.
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 two-phase flows in a typical 3D model pore network of the chamber-and-throat 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 two-phase flows in porous media in mixtures of connected and disconnected flows (connected-oil pathway flow, dynamics of small and large ganglia, drop traffic flow).
Table 3 Reduced pressure gradient values, x, for different flow conditions, (Ca, r), and viscosity ratios, k.Systematic mutations have been revealed within the flow with changing flow conditions.Some of them surface-up 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 connected-disconnected flows.Deeper within the disconnected-oil flow, strong/ intense mutations between Ganglion Dynamics and Drop-Traffic 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.This characteristic behavior is attributed to the interstitial balance between capillarity and bulk viscosity.Flow of the non-wetting phase in disconnected form is significant and for certain flow conditions the dominant flow mode.
In general, sustainability of two-phase 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 (k ≠ 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 non-linearity 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 non-wetting phase flow.
Concluding, the present work indicates the potential of the DeProF true-to-mechanism theoretical framework in revealing the effects of underlying flow mechanisms in immiscible steady-state two-phase 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.
all unit cells) surface area of oil/water interfaces, advancing and receding P all unit cells) surface area of o/w interfaces (mobilized and stranded) of any iall unit cells) surface area of the mobile oil/water interface of any mobilized ganglion S all unit cells) surface area of the oil/water interface of any i-class stranded ganglion S (30)ÃsMG ow;iMean (across all unit cells) surface area of the stranded oil/water interface of any i-class mobilized ganglion S Size of the ensemble of physically admissible solutions M (16) W Rate of mechanical energy dissipation within the process M (56) W 1F Mechanical power dissipation for an equivalent 1-phase flow M (56) W C Mechanical power dissipated in CPF M (57)

Figure 1 (
Figure 1 (a) Schematic representation of the actual flow (left sketch) and its theoretical decomposition into prototype flows: Connected-oil Pathway Flow (CPF) and Disconnected-Oil 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 unit-cells.The thick dashed line separates the Ganglion Dynamics (GD) cells domain and the Drop-Traffic Flow (DTF) cells domain.

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 two-phase flow visits a canonical ensemble of physically admissible flow configurations S 0 w ; b 0 ; v 0 depicted by the cloud of small, red balls.Here, Ca = 1.19 Â 10 À6 , k = 1.45 and flow rate ratio values, r, span 2 orders of magnitude.A unique triplet of values (S w , b, v), 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.

Figure 4
Figure 4Typical jik-class 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.

Figure 5
Figure 5 Input/output block diagram in the DeProF model implementation.The variables describing the structure of the interstitial flow are presented in Section 3.
(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 ℓ ffiffi ffi 2 p and surface area equal to ffiffi ffi 3 p ℓ 2 =4.When the cubic lattice elements are stacked, these equilateral triangles stick side-by-side and exactly fill the frontal area.An impression of a 3Â2Â1 stacking of unit-lattice 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 unit-cells aligned to the macroscopic flow, are virtually confined within a normal equilateral twisting prism.The prism sides are equal to ℓ ffiffi ffi 2 p

Figure 6 (
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 checker-board pattern and exactly fill any frontal area perpendicular to the macroscopic flow.

Figure 7 A
Figure 7 A close-up 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.
Figure 9Mean 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, m G i ; i ¼ 1; 3; 5; 7; 9.

Figure 8 (
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 end-gate 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.J represents the idealized radii of curvature of the o/w interfaces for advancing and receding contact angles.
as for o/w interfaces in the DTF domain, A owMD , A owSD , A owD .
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 disconnected-oil flow domain, U 0 o;DOF , as the former comprise an inseparable mixture within DOF.They are given by for different values of the capillary number, the flow rate ratio and the viscosity ratio, are presented in Section 4.

Figure 10 Expected
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 k ¼ mo =m w , Arrows indicate increasing Ca values: !10 À8 Ca 10 À4 and ---> 4 Â 10 À7 Ca 10 À4 .

Figure 11 Diagrams
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, ðn 0 1 ; n 0 2 ; z 0 Þ.The values of flow conditions, (Ca, r), and average flow arrangement variables,(S w , b, v), are provided in the upper two rows of each tablet.
Figure 12 Expected, ensemble average values of the flow arrangement variables for different values of the capillary number, Ca, and viscosity ratio, k.(a) Saturation, S w ; (b) volume fraction of connected pathway flow (CPF), b; (c) volume fraction of ganglion dynamics flow (GD) within the disconnected oil flow (DOF), v. Arrows indicate increasing Ca values: !10 À8 Ca 10 À4 and ---> 4 Â 10 À7 Ca 10 À4 .

4. 2
Typical, physically admissible flow configurations Indicative physically admissible flow configurations described by the values of the FAV, for different flow conditions, are presented in the diagrams of Figure 11.'Moderate flow arrangements' are characterized by physically admissible flow arrangement values ðS 0 w ; b 0 ; v 0 Þ close to the ensemble average, (S w , b, v) and are presented in diagrams (a1)-(c1).'Extreme flow arrangements' are characterized by physically admissible flow arrangement values ðS 0 w ; b 0 ; v 0 Þ pertaining to the extreme corners of the ensemble envelope and are presented in diagrams (a2)-(c2).

4. 3
Flow arrangement variables (steady-state average)Expected values of the FAV, i.e. the saturation, S w , the volume fraction of CPF, b, the volume fraction of GD within the DOF, v, for different values of the capillary number, Ca, and the flow rate ratio, r, for different viscosity ratio systems, k, are presented in the diagrams of Figure12.As the viscosity ratio is increased, the onset and build-up of connected oil pathways (b) are observable from lower Ca values.

Figure 14 Expected
Figure 14 Expected, ensemble average values of the mean ganglion sizes for different values of the capillary number, Ca, and viscosity ratio, k. a) Mean size of all ganglia (mobilized þ stranded), V G m ; b) mean size of mobilized ganglia, V mG m .Arrows indicate increasing Ca values: !10 À8 Ca 10 À4 and ---> 4 Â 10 À7 Ca 10 À4 .

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, k ¼ mo =m w 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, 10-fold 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, Ãm ow , whether these are ganglia, ÃmG ow , or droplets, ÃmD ow , so that Ãow ¼ Ãm ow þ Ãs ow and Ãm ow ¼ ÃmG ow þ ÃmD ow .The significant drop in the overall o/w interface area density, A ow , observed in Figure 15d at intense two-phase flow conditions (high Ca, high r values), is due to the increase of the volume fractions of CPF, see Figure 12b .

Figure 16 Expected
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, k ¼ mo =m w .Arrows indicate increasing Ca values: !10 À8 Ca 10 À4 and ---> 4 Â 10 À7 Ca 10 À4 .

4. 7
Dynamic variables describing the interstitial transport properties of the process Ensemble average values of the fractional flow of oil through CPF, h o,C , GD, h o,G , and DTF, h 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.
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, k ¼ mo =m w .a) Reduced superficial velocity of o/w surface, U ow,DOF , b) Fraction of o/w surface transport through ganglion dynamics, j ow.G , c) Fraction of o/w surface transport through drop traffic flow, j ow.D .Arrows indicate increasing Ca values: !10 À8 Ca 10 À4 and ---> 4 Â 10 À7 Ca 10 À4 .
while the vertical, dashed line in red indicates the corresponding flow rate ratio, r Ã ∞ ¼ 1= ffiffi ffi k p

Figure 19 Effect
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, k ¼ mo =m w .Arrows indicate increasing Ca values: !10 À8 Ca 10 À4 and ---> 4 Â 10 À7 Ca 10 À4 .

Figure 20 The
Figure 20 The set of critical flow conditions, r * i ; Ca i n o , detected by the DeProF simulations, for which the energy efficiency takes maximum values (ref.Fig. 19).

Table 1
Occurrence probabilities, f j , and respective reduced size diameters for the 5 classes of chambers, D Cj , and throats, D Tj , of the 3D pore network used in the DeProF simulations.All dimensions are normalized to the lattice constant, ℓ = 1221 mm.

Table 2
Volume fraction of the connected pathway flow (CPF) domain, b, against fractional flow of oil through it, h o,C , for various flow conditions, (Ca, r), and viscosity ratios k.discrepancies are pronounced in the low Ca / low r regimes and they to fade-away 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.