© M.S. Valavanides, published by IFP Energies nouvelles, 2018

Notation

∼: A tilde embellishment indicates a dimensional variable

No tilde embellishment denotes a dimensionless variable

Φ^′: A prime indicates a physically admissible value of variable Φ (see Sect. 1.5)

In the following lists the far right columns indicate:

Scale: M (macroscopic), S (statistical average), P (pore scale) (see Sect. 1.2, end)

Definition: corresponding equation (number) or paragraph number

Latin characters

Notation	Interpretation of variable	Scale	Definition
${\displaystyle \tilde{A}}$	Porous medium cross-sectional area of REV	M	(2)
${A^{\prime }}_{ow}$	Total o/w surface area	S	(43)
${A^{\prime }}_{owD}$	Total o/w surface area in the DTF domain	S	(41)
${A^{\prime }}_{owDOF}$	Total o/w surface area in the DOF domain	S	(42)
${A^{\prime }}_{owG}$	Total o/w surface area of ganglia	S	(36)
${A^{\prime }}_{owMD}$	Oil/water surface area of mobilized drops, DTF domain	S	(39)
${A^{\prime }}_{owMG}$	Oil/water surface area of mobilized ganglia	S	(33)
${A^{\prime }}_{owMMG}$	Mobilized o/w surface area of mobilized ganglia	S	(31)
${A^{\prime }}_{owSD}$	Oil/water surface area of stranded drops, DTF domain	S	(40)
${A^{\prime }}_{owSG}$	Oil/water surface area of stranded ganglia	S	(34)
${A^{\prime }}_{owSMG}$	Stranded o/w surface area of mobilized ganglia	S	(32)
$A_{ow}^{adv}$, $A_{ow}^{rec}$	Mean (across all unit cells) surface area of oil/water interfaces, advancing and receding	P	(27)
${{\displaystyle \tilde{A}}}_{ow,i}^{MG}$	Mean (across all unit cells) surface area of o/w interfaces (mobilized and stranded) of any i-class mobilized ganglion	S	(29)
${{\displaystyle \tilde{A}}}_{ow}^{mMG}$	Mean (across all unit cells) surface area of the mobile oil/water interface of any mobilized ganglion	S	(27)
${{\displaystyle \tilde{A}}}_{ow,i}^{SG}$	Mean (across all unit cells) surface area of the oil/water interface of any i-class stranded ganglion	S	(30)
${{\displaystyle \tilde{A}}}_{ow,i}^{sMG}$	Mean (across all unit cells) surface area of the stranded oil/water interface of any i-class mobilized ganglion	S	(28)
Ca	Capillary number	M	(1)
${C^{\prime }}_i^{C,G}$	Fraction of class-i core (‘C’) ganglion cells − conducting	S	(20)
${C^{\prime }}_i^{X,G}$	Fraction of class-i non-conducting (‘X’) ganglion cells	S	(21)
${C^{\prime }}_i^{E,G}$	Fraction of class-i end (‘E’) ganglion unit cells − conducting	S	(22)
DCj	Diameter of j-class chamber	P	2.1
DTi	Diameter of i-class throat	P	2.1
fEU	Energy efficiency of the process	M	(58)
fj	Occurrence probability of any i-class network element	P	2.1
$f_{OF}^{\text{'}}$	Coefficient of oil fragmentation	S	(48)
${g^{\prime }}_{a,i}^b$, ${g^{\prime }}_a^D$	Conductance of the various unit cells in the GD and DTF domains	P	1.3
Imax	Maximum ganglion size class	P	1.4
${\displaystyle \tilde{k}}$	Absolute permeability of the pore network	M	(5)
ℓ	Network lattice constant	P	2.1
mD	Mean mobilization probability of drops	S	2.3
$m_i^G$	Mean mobilization probability of i-class ganglia	S	2.3
$n_i^G$	Number density of i-class ganglia	S	(15)
Nd	Network dimensionality	P	2.1
$N_i^G$	Number of pore unit cells occupied by each i-class ganglion	P	(18)
${{\tilde{\mathit{N}}}_{\mathit{D}}}^{\prime }$	Average number of oil drops p.u.v. of DTF domain	S	(38)
${N^{\prime }}_{crit}^D$	Critical number density of minimum size oil droplets	P	(46)
$P_{a,i}^b$, $P_a^D$	Occurrence probabilities of unit cells in the GD and DTF domains	S	(13)
${{\displaystyle \tilde{q}}}_o$, ${{\displaystyle \tilde{q}}}_w$	Macroscopic flow rates of oil and water	M	(2)
${{\displaystyle \tilde{q}}}_{o,C}$	Mean flow rate of oil through CPF	M	(49)
${{\displaystyle \tilde{q}}}_{uc,i}^G$	Flow rate within any conducting i-class ganglion unit cell	P	(19)
r	Oil/water flow rate ratio	M	(2)
Sw	Mean water saturation (entire network)	M	1.2
So,D	Mean oil saturation in the DTF domain	S	(9)
So,G	Mean oil saturation in the GD domain	S	(9)
$u_i^G$	Mean ganglion velocity of i-class ganglia	S	(12)
${U^{\prime }}_{o,D}$	Superficial velocity of oil in the DTF domain	M	(50)
${U^{\prime }}_{o,DOF}$	Superficial velocity of oil in the DOF domain	S	(12)
${U^{\prime }}_{o,G}$	Superficial velocity of oil in the GD domain	S	(50)
${U^{\prime }}_{o,C}$	Superficial velocity of oil in the CPF	M	(6)
${{\displaystyle \tilde{U}}}_o$,${{\displaystyle \tilde{U}}}_w$	Superficial velocity of oil and water	M	(2)
${U^{\prime }}_{ow,G}$	Superficial velocity of the o/w interface in the GD domain	S	(52)
${U^{\prime }}_{ow,mD}$	Superficial velocity of the o/w surface in the DTF domain	S	(53)
${U^{\prime }}_{ow,DOF}$	Superficial velocity of the o/w interface in the DOF domain	M	(54)
Uw	Superficial velocity of water	M	(2)
Uw,DOF	Superficial velocity of water in the DOF domain	M	(7)
${\displaystyle \tilde{V}}$	Pore network reference volume	M	1.2
VD	Average (mean) oil drop volume	S	(37)
VD,1	Minimum size of an oil-droplet	P	(44)
$V_i^G$	Volume of i-class ganglia pore unit cells	P	(24)
$V_m^G$	Mean size of all ganglia (mobilized and stranded)	S	(24)
$V_m^{MG}$	Mean size of the mobilized ganglia	S	(25)
${V^{\prime }}_m^{SG}$	Mean size of the stranded ganglia	S	(26)
$V_{{\mathrm{\Omega }}_{PAS}}$	Size of the ensemble of physically admissible solutions	M	(16)
W	Rate of mechanical energy dissipation within the process	M	(56)
${{\displaystyle \tilde{W}}}^{1\mathrm{\Phi }}$	Mechanical power dissipation for an equivalent 1-phase flow	M	(56)
WC	Mechanical power dissipated in CPF	M	(57)
WD	Mechanical power dissipated in DTF	M	(57)
WG	Mechanical power dissipated in GD	M	(57)
x	Macroscopic pressure gradient	M	(5)
xpm	Set composed of pore network parameters	P	1.1
${\displaystyle \tilde{z}}$	Direction of macroscopic flow	M	(5)

Greek characters

Notation	Interpretation of variable	Scale	Definition
β	Porous medium volume fraction occupied by the connected oil	M	(8)
ζ	Parameter measuring the sharpness of the ganglion size distribution	S	(15)
$\eta {\text{'}}_{o,C}$	Fractional flow through connected pathways	M	(49)
$\eta {\text{'}}_{o,D}$	Fractional flow through drop traffic flow	M	(49)
$\eta {\text{'}}_{o,G}$	Fractional flow through ganglion dynamics	M	(49)
θA, θR	Advancing and receding contact angles of the o/w interface	P	1.1
κ	Oil/water viscosity ratio	M	(3)
${\stackrel{\sim }{\mathit{\gamma }}}_{\mathit{o}\mathit{w}}$	Oil/water interfacial tension	P	(1)
${\stackrel{\sim }{\mathit{\mu }}}_{\mathit{w}}$	Viscosity of water	M	(3)
${\stackrel{\sim }{\mathit{\mu }}}_{\mathit{o}}$	Viscosity of oil	M	(3)
$\mathit{\xi }{\text{'}}_{ow,D}$	Fractional transport of the o/w interfaces through DTF	S	(55)
$\mathit{\xi }{\text{'}}_{ow,\mathit{G}}$	Fractional transport of the o/w interfaces through GD	S	(55)
σ3D	Mean coordination number of pore network	P	2.1
χG	Tortuosity of ganglia in the pore network	P	2.1
ω	Fraction of all the ganglion cells over all the DOF region cells	M	(9)

Abbreviations

CPF	Connected Pathway Flow
DeProF	Decomposition in Prototype Flows
DOF	Disconnected Oil Flow
DTF	Drop Traffic Flow
FAV	Flow Arrangement Variable(s)
GD	Ganglion Dynamics
GCD	Ganglion Cells Domain
o/w	oil/water
p.m.	porous medium
p.u.v.	per unit volume
p.u.v.p.m.	per unit volume of porous medium
REV	rRepresentative Elementary Volume
u.c.	unit cell

Introduction

Fluid flow in porous media, i.e. porous rocks, fracture networks and granular media, is subject to considerable interdisciplinary research activity in physics, earth sciences and engineering. Examples of natural and engineered processes include hydrocarbon recovery, carbon dioxide 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, non-wetting, 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 and Payatakes, 1996 and 2000; 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 and Payatakes, 1995 and 1999; Krummel et al., 2013; Armstrong et al., 2016). These studies show that during simultaneous flow of oil and water, the disconnected phase (oil) contributes significantly to the total flow and, in certain cases even exclusively, e.g. emulsion flow (Cobos et al., 2009; Guillen et al., 2012). Furthermore, the flow rate against pressure gradient relation is found to be strongly 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, non-wetting 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 two-phase flows in pore networks, Valavanides and Payatakes (2000 and 2001; also Valavanides, 1998 and 2012) developed the DeProF mechanistic model. The model is based on the concept of decomposition in prototype flows (see next section). It accounts the 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₂ 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.

1 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.

1.1 The independent variables for steady-state two-phase 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, ${{\displaystyle \tilde{U}}}_o$ and ${{\displaystyle \tilde{U}}}_w$, conventionally defined as oil and water flow rates, ${{\displaystyle \tilde{q}}}_o$ and ${{\displaystyle \tilde{q}}}_w$, divided by the p.m. cross-sectional area, ${\displaystyle \tilde{A}}$. The set of superficial velocities may be appropriately reduced and replaced by a set of dimensionless variables, comprising the capillary number, $$Ca=\frac{{{\displaystyle \tilde{\mu }}}_w{{\displaystyle \tilde{U}}}_w}{{{\displaystyle \tilde{\gamma }}}_{ow}}$$(1) and the oil/water flow rate ratio, $$r=\frac{{{\displaystyle \tilde{q}}}_o}{{{\displaystyle \tilde{q}}}_w}=\frac{{{\displaystyle \tilde{U}}}_o{\displaystyle \tilde{A}}}{{{\displaystyle \tilde{U}}}_w{\displaystyle \tilde{A}}}=\frac{{{\displaystyle \tilde{U}}}_o}{{{\displaystyle \tilde{U}}}_w}$$(2)where ${{\displaystyle \tilde{\mu }}}_w$ is the viscosity of water and ${{\displaystyle \tilde{\gamma }}}_{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., ${\displaystyle \tilde{k}}$_;

• the oil/water viscosity ratio;

$$\mathit{\kappa }={{\displaystyle \tilde{\mathit{\mu }}}}_o/{{\displaystyle \tilde{\mathit{\mu }}}}_w$$(3)

• the dynamic advancing and receding contact angles of the o/w interface on the pore walls, θ_A and θ_R, and;

• a (sub)set, x_pm, comprising all the dimensionless geometrical and topological parameters of the pore network that affect the flow (e.g. porosity, genus, coordination number, normalized chamber and throat size distributions, chamber-to-throat size correlation factors, etc).

1.2 The concept of Decomposition into Prototype Flows (DeProF)

In the general case of steady-state, one-dimensional, two-phase 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, ${\displaystyle \tilde{z}}$, with constant flow rates, ${{\displaystyle \tilde{q}}}_o$ and ${{\displaystyle \tilde{q}}}_w$, 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).

All types of flow have been observed experimentally in model porous media (Avraam and Payatakes, 1995, 1999; Tsakiroglou et al., 2007; Tallakstad et al., 2009; Krummel et al., 2013; Armstrong et al., 2016) as well as in real porous media (Van de Merwe and Nicol, 2009; Georgiadis et al., 2013; Oughanem et al., 2015; Rücker et al., 2015).

In the DeProF model, each prototype flow (CPF, GD and DTF) has the essential characteristics of the corresponding flow pattern in suitably idealized form and so, the pore scale mechanisms are incorporated in the prototype flows. In the Connected Pathway Flow (CPF) region the oil retains its connectivity and flows with virtually one-phase flow. The p.m. volume fraction occupied by the connected oil is denoted by β. 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 − β). 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 ${\displaystyle \tilde{V}}$. The size (volume), ${\displaystyle \tilde{V}}$, 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 ω, and is called the GD domain fraction. The DTF domain fraction in the DOF region equals (1 − ω).

S_w, β and ω are called Flow Arrangement Variables (FAV) because they provide a coarse indication of the structure (or ‘immiscible mixing’ of the constituent prototype flows) of the flow pattern. One of the objectives of the DeProF model algorithm is to determine the values of S_w, β and ω that conform with the externally imposed flow conditions (Ca, r).

The flow analysis is carried-out on two length scales, the macroscopic scale (10¹² 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.

1.3 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 $$x=x(Ca,r;\kappa ,{\theta }_A^0,{\theta }_R^0,x_{pm})$$(4) where, $$x=-\frac{\partial {\displaystyle \tilde{p}}}{\partial {\displaystyle \tilde{z}}}\frac{{\displaystyle \tilde{k}}}{{{\displaystyle \tilde{\gamma }}}_{ow}Ca}$$(5)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 ${{\displaystyle \tilde{U}}}_w$ [the 2nd component of the product in Equation (5)].

Referring to Figure 1, the externally imposed flow rates of oil, ${{\displaystyle \tilde{q}}}_o$, and water, ${{\displaystyle \tilde{q}}}_w$, are partitioned in the prototype flows. In particular, ${{\displaystyle \tilde{q}}}_{o,C}$, ${{\displaystyle \tilde{q}}}_{o,DOF}$ are the mean oil flow rates in CPF and DOF and ${{\displaystyle \tilde{q}}}_{w,DOF}$ the mean water flow rate in DOF (as CPF is oil-saturated, ${{\displaystyle \tilde{q}}}_{w,DOF}={{\displaystyle \tilde{q}}}_w$). 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}={{\displaystyle \tilde{q}}}_{o,C}/{{\displaystyle \tilde{q}}}_o$, and the reduced oil and water mean superficial velocities defined by $U_{m,DOF}={{\displaystyle \tilde{q}}}_{m,DOF}/{{\displaystyle \tilde{q}}}_m$, where m = o, w.

The flow analysis produces the following system of equations (Valavanides and Payatakes, 2000; 2001) $$U_{{}^{o,C}}=\frac{1}{\kappa r}x$$(6) $$\left(1-\beta \right)U_{{}^{w,DOF}}=1$$(7) $$\beta U_{{}^{o,C}}+\left(1-\beta \right)U_{{}^{o,DOF}}=1$$(8) $$\beta +\left(1-\beta \right)\left[\omega S_{{}^{o,G}}+\left(1-\omega \right)S_{{}^{o,D}}\right]=1-S_{{}^w}$$(9) $${\displaystyle \sum _{i=1}^{I_{\max }}}{n}_i^G{N}_i^G=2N_d$$(10) $$\frac{1}{2}{\displaystyle \sum _{i=1}^{I_{\max }}}i{n}_i^G=S_{{}^{o,G}}$$(11) $${\displaystyle \sum _{i=1}^{I_{\max }}}i{n}_i^G{m}_i^G{u}_i^G=2r{U}_{{}^{o,DOF}}$$(12) $$\omega {\displaystyle \sum _{a,b,i}}{P}_{a,i}^b\frac{\left(1+r\right)-x\beta /\kappa -x{g}_{a,i}^b\left(x\right)\left(1-\beta \right)/k}{\left(1+r-x\beta /\kappa \right)+x{g}_{a,i}^b\left(x\right)\left(1-\beta \right)/k}+\left(1-\omega \right){\displaystyle \sum _a^{\mathrm{}}}{P}_a^D\frac{\left(1+r\right)-x\beta /\kappa -x{g}_a^D\left(x\right)\left(1-\beta \right)/k}{\left(1+r-x\beta /\kappa \right)+x{g}_a^D\left(x\right)\left(1-\beta \right)/k}=0$$(13) where _{$n_i^G,i=\mathrm{1,...},I_{\max }$} is the reduced number density of the i-class ganglia and I_max is the upper truncation ganglion size class. _{$N_i^G$} is the number of unit cells occupied by any i-class ganglion. _{$u_i^G(={{\displaystyle \tilde{u}}}_i^G/{{\displaystyle \tilde{U}}}_w)$} is the reduced mean ganglion velocity i-class ganglia. _{$m_i^G$} is the mean mobilization probability of the i-class ganglia. _{$P_{a,i}^b$} and $P_a^D$ are the occurrence probabilities of the various unit cells in the GD and DTF domains respectively. _{$g_{a,i}^b$} and $g_a^D$ are the reduced conductance of the various unit cells in the GD and DTF domains respectively.

Note that m^D, $u_i^G$, $m_i^G$, $g_{a,i}^b$ and $g_a^D$ depend on x, as described by semi-analytical expressions of momentum balance derived over the pore scale. $P_{a,i}^b$, $P_a^D$, $g_{a,i}^b$ and $g_a^D$ 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_a^D$ and $g_a^D$ 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, ${\displaystyle \tilde{V}}$; 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 freedom by allowing the process to attain various ganglion size distributions and deliver a variety of physically admissible solutions − see Section 1.5. Equation (13) results from application of Effective Medium Theory (Kirkpatrick, 1973) to the equivalent 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 (I_max + 6) dependent variables, namely x, U_o,CPF, U_o,DOF, U_W,DOF, S_o,G, S_o,D and $\left\{n_i^G,i=1,\mathrm{...},I_{\max }\right\}$, along with the three FAV, namely S_w,  β,  ω. The system is closed by imposing an appropriate type of (sharply decreasing) distribution function for the ganglion volumes, Equation (15), that is based on modeling assumptions presented in the following paragraph. This way the system reduces to a single algebraic equation, as Equations (6)–(12) are all 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_{a,i}^b$, $P_a^D$, $g_{a,i}^b$ and $g_a^D$ on the given pore network geometry. This may be accomplished either by furnishing appropriate analytical expressions (if possible) for $P_{a,i}^b$, $P_a^D$, $g_{a,i}^b$ and $g_a^D$ 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.

1.4 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 β region of the pore network and remain aligned to the macroscopic flow. Of a total number of M unit-cells in a reference volume, βM unit-cells are occupied by connected oil (the ensemble of Connected-Oil Pathways) and allow the rest (1 − β) 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 − β)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 − β) ω of the pore network, leaving to the DTF cells a complementary volume fraction of (1 − β) (1 − ω). The aforementioned partitioning is described by Equations (7)–(9). The DTF unit cells are all considered to be indistinguishable − refer to Section 2.2 and Figure 7.

Ganglion 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 $${{\displaystyle \tilde{n}}}_i^G=\frac{\text{number}\text{of}\mathit{\text{i}}-\text{class}\text{ganglia}}{\text{GCD}\text{volume}}=n_i^G\frac{1}{2N_d{{\displaystyle \tilde{V}}}_{uc}}$$(14) where N_d is the network dimensionality (=3 for 3D networks), ${{\displaystyle \tilde{V}}}_{uc}$ is the average volume of unit cells and $n_i^G$ 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 $$n_i^G=\{\begin{array}{ccc}\hfill n_i^G\hfill & \hfill ,\hfill & \hfill i=1,2\hfill \\ \hfill n_i^G{\zeta }^{\left(i-2\right)}\hfill & \hfill ,\hfill & \hfill 3\le i\le I_{\max }\hfill \\ \hfill 0\hfill & \hfill ,\hfill & \hfill i>I_{\max }\hfill \end{array}$$(15)where n₁,  n₂,  I_max and ζ are parametric values that define the particular physically admissible flow configuration. I_max is the maximum attainable size of a ganglion and 0 < ζ < 1 provides a measure of the sharpness of the ganglion size distribution, see Figure 2.

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., 20142; Youssef et al., 2014; Georgiadis et al., 2013; Datta et al., 2014; Oughanem et al., 2015; Armstrong et al., 2016).

Indicative DeProF model predictions of physically admissible ganglion size distributions, $(n_2^G,n_2^G,\zeta )$, are presented in the diagrams in Section 3.1, along with the corresponding values (tabulated) of the FAV (S_w, β, ω). The expected, average values of the flow arrangement variables, (S_w, β, ω), for different values of the capillary number, the flow rate ratio and the viscosity ratio, are presented in the diagrams in Section 3.2.

1.5 Physically Admissible Flow Configurations

A core feature of the DeProF model is the detection of all flow configurations that are physically admissible under the imposed macroscopic flow conditions. The actual flow at any REV of the p.m. (extending at an order of 10⁹ pores along the main flow direction), ‘wanders’ within the domain of physically admissible flow configurations, ‘visiting’ any one with equal probability or frequency. This equiprobability is a fundamental assumption of the DeProF model, related to the ergodic nature of the process (see ‘2.3. Discussion on Assumption of the Process Ergodicity’ in Valavanides and Daras, (2016)

To detect these flow configurations, the domain {S_w, β, ω} of all possible flow arrangement values is partitioned using sufficiently fine steps to obtain a 3D grid (Fig. 3). Then, for each combination of the flow arrangement values (S_w, β, ω), the DeProF algorithm is asked to detect (if) a solution of the DeProF Equations (6)–(13) and (15) (exists). In such case the selected triple of {S_w, β, ω} values, is allowed as a physically admissible solution (PAS) and is denoted as $\left\{{S^{\prime }}_w,{\beta }^{\prime },{\omega }^{\prime }\right\}$. Otherwise, the triple is rejected. At the end of the discretization and scanning process, the sub-domain of the {S_w, β, ω} space, denoted Ω_PAS, is detected as the set of physically admissible solutions. This set is depicted by the cloud of small red balls in the diagrams of Figure 3.

The measure (“volume”) of Ω_PAS is given by $$V_{{\mathrm{\Omega }}_{PAS}}\left(Ca,r\right)={\displaystyle \underset{{\mathrm{\Omega }}_{PAS}\left(Ca,r\right)}{\iiint }}d{S}_{w}d\beta d\omega$$(16) where Ω_PAS(Ca, r) stands for integration carried over the physically admissible ranges in {S_w, β, ω} for the imposed pair of (Ca, r) values. The volume of Ω_PAS is a measure of the degrees of freedom of the process and the variety of microstates it can attain; it is also related to the rate of entropy production (configurational entropy) within the process (Valavanides and Daras, 2016). The volume of the ensemble of interstitial flows is a measure of the internal degrees of freedom of the process; its size and shape depend on the externally imposed macroscopic flow conditions. There are flow conditions for which there are many physically admissible interstitial flow arrangements and there are flow conditions for which the macroscopic flow is confined within a tight ensemble: note the change in position, size and shape of the red cloud with changing flow conditions, Figure 3.

By assuming that each physically admissible solution is visited with the same probability, or frequency (assumption of ergodicity), and averaging over their domain, Ω_PAS, a unique solution for the macroscopic flow is obtained. For any quantity, Φ ', the corresponding expected mean macroscopic flow quantity, Φ is defined as $$\mathrm{\Phi }\left(Ca,r\right)=\langle \mathrm{\Phi }\text{'}\rangle =\frac{{\displaystyle \underset{{\mathrm{\Omega }}_{PAS}\left(Ca,r\right)}{\iiint }}\mathrm{\Phi }\text{'}\left(Ca,r\right)d{S}_{w}d\beta d\omega }{{\displaystyle \underset{{\mathrm{\Omega }}_{PAS}\left(Ca,r\right)}{\iiint }}d{S}_{w}d\beta d\omega }=\frac{1}{V_{{\mathrm{\Omega }}_{PAS}}\left(Ca,r\right)}{\displaystyle \underset{{\mathrm{\Omega }}_{PAS}\left(Ca,r\right)}{\iiint }}\mathrm{\Phi }\text{'}\left(Ca,r\right)d{S}_{w}d\beta d\omega$$(17)

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, β, ω) can be obtained by averaging over the domain of physically admissible solutions, using Equation (17), with Φ replaced by S_w, β and ω. These values define the ensemble average flow configuration of the macroscopic flow. They are indicated by the large black balls in the diagrams of Figure 4.

In computational practice, the size of the domain of physically admissible solutions, $V_{{\mathrm{\Omega }}_{PAS}}$, as well as any expected mean macroscopic flow quantity, Φ, can be determined by any coarse graining (discretization) procedure used to detect the ensemble of physically admissible solutions and delineate its envelope. The phase space {S_w × β × ω} pertaining to the FAV, is partitioned in corresponding N_S × N_β × N_ω phase space voxels. The finer the graining, the smaller becomes the size of these voxels. Nevertheless, the smallest discretization size of voxels has a limit; it cannot become smaller than a minimum, critical size in each prototype flow. Suppose that we detect two solutions ${({S^{\prime }}_w,{\beta }^{\prime },{\omega }^{\prime })}_i$ and ${({S^{\prime }}_w,{\beta }^{\prime },{\omega }^{\prime })}_{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 re-discretization 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, $({{\displaystyle \tilde{U}}}_o,{{\displaystyle \tilde{U}}}_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, submitted for review). Additional comments are provided in Section 4.9 ‘Discussion of results’ in the present work.

2 Pore Network

Implementation of the DeProF analysis and methodology presented in the previous paragraphs is directly related to the structure of the pore network and the associated modeling provisions pertaining to the prototype flows and, especially, to the disconnected fluidic elements. The following paragraphs describe the essential characteristics of the typical model pore network used in the present simulations.

2.1 Pore Network Geometry

The simulations are performed on the 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 and Payatakes (1995, 1999).

The 3D cubic lattice pore network is of the “ball-and-stick” 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 σ_3D = 6. For the particular network dimensions, the volume porosity is ϵ = 0.0464 and its absolute permeability is ${\displaystyle \tilde{k}}=8.965\mu {\text{m}}^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 $\ell \sqrt{2}/2$ and volume ${\ell }^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(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 $\ell \sqrt{2}$ and surface area equal to $\sqrt{3}{\ell }^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 $\ell \sqrt{2}$ and its frontal area is equal to $\sqrt{3}{\ell }^2/4$.

An important parameter associated with the ‘zigzag' motion of the ganglion is the tortuosity of its spine, χ_G. This is also the ratio of the length of its actual migration path over the effective migration length in the macroscopic flow direction. This parameter is essentially one of the p.m. geometrical and topologic characteristics and will be considered to be equal for all ganglion size classes. For the case under consideration, ${\chi }_G={\sin }^{-1}\arctan \left(\sqrt{2}/2\right)\approx 1.732$.

2.2 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_i^G$of unit cells occupied by each i-class ganglia are given by $$N_i^G=i\left(2N_d-1\right)+1$$(18) where N_d is the dimensionality of the network lattice. For the 3D lattice used in the simulations N_d = 3 (but for the 2D depiction of the lattice in Figure 3, N_d = 2)

Unit cells in which ganglion menisci exist are called gate unit cells. In the present work, we assume that all ganglia attain a ‘wormlike’ alignment − as much as the 3D zigzag pore conduits allow − to the macroscopic flow direction. Among all gate unit cells of a mobilized ganglion, there are only two gate unit cells in which rheons take place (one cell for the ‘xeron’ and one cell for the ‘hygron’ (Payatakes, 1982). Such cells are denoted with the index ‘E, G’ (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, ${{\displaystyle \tilde{q}}}_{uc,i}^G$, 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, ${{\displaystyle \tilde{u}}}_i^G$, i.e., $${{\displaystyle \tilde{q}}}_{uc,i}^G={{\displaystyle \tilde{u}}}_i^G{\ell }^2{\chi }_G$$(19) where χ_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_i^G$ 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, $${C^{\prime }}_i^{C,G}=m_i^G(i-1)/N_i^G$$(20) Non-conducting ganglion cells (note here, that this fraction takes into account all the cells of stranded ganglia), $${C^{\prime }}_i^{X,G}=\left[i\left(2N_d-1\right)+1-m_i^G\left(i+1\right)\right]/N_i^G$$(21)End ganglion cells − conducting, through which hygrons take place, $${C^{\prime }}_i^{E,G}=2m_i^G/N_i^G$$(22)Now, as expected, $${C^{\prime }}_i^{C,G}+{C^{\prime }}_i^{X,G}+{C^{\prime }}_i^{E,G}=1$$(23)

Drop Traffic Flow unit cells

The unit cells in the DTF domain contain water and oil droplets. Water may be connected or not depending on the formation of o/w menisci in contact with the pore walls. This is a modeling provision, associated − in general − to the wetting characteristics of the system, i.e. the combination of the dynamic contact angles, the size of the droplets and the pore wall 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 droplet-sizes. 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.

2.3 Mobilization of drops and ganglia

Oil droplets and ganglia have a tendency to get stranded, due to the capillary pressure difference that can be maintained across their o/w interfaces. This is attributed to the hysteresis between the advancing/receding contact angles that, in combination with the local pore geometry configuration, result in the hysteresis of the mean curvatures between the advancing and receding menisci. The mobilization probability of a ganglion or a droplet depends on the microscopic pressure difference − locally induced − along every pair of o/w interfaces and the net capillary pressure differences over the respective o/w interfaces (Payatakes, 1982).

In the DeProF model, stranding of ganglia and droplets is examined over particular stranding configurations, i.e. layouts within the unit cells they reside in, for which they get stranded. For these stranding configurations within the unit cells the capillary pressure differences become large and stable. Stability means that if any oil blob (drop or ganglion) is slightly perturbed, the new capillary pressure differences become larger and force the blob to return into its original position, keeping it stranded, until − eventually − another, larger pressure difference (over) perturbs the blob and destabilizes it; then it migrates downstream in consecutive episodic rheons, triggering Haines jumps and propagating pressure pulses (Berg et al., 2013; Rücker et al., 2015). For each type of blob and possible configuration, a respective mobilization probability is assigned to a solitary blob, i.e. a single blob in the pore network. Then, to account the crowding effect of the motion of neighboring blobs (whether drops or ganglia) induced by their episodic pore scale events, the locally induced macroscopic pressure gradient is perturbed by imposing a 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.

3 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.

3.1 Mean ganglion sizes

The expected values of the mean ganglion size are computed after a physically admissible solution, specifically, a physically admissible ganglion size distribution, ${n^{\prime }}_i^G,i=\mathrm{1,...},{I^{\prime }}_{\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^{\prime }}_m^G$, is estimated as the ratio of the entire oil ganglia mass over the entire ganglia population, as $${V^{\prime }}_m^G=\frac{{\displaystyle \sum _{i=1}^{I_{\max }}}{n^{\prime }}_i^G{V}_i^G}{{\displaystyle \sum _{i=1}^{I_{\max }}}{n^{\prime }}_i^G}$$(24) where $V_i^G$ 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^{\prime }}_m^{MG}$, is estimated as the ratio of the entire oil mass of mobilized ganglia over the population of mobilized ganglia, as $${V^{\prime }}_m^{MG}=\frac{{\displaystyle \sum _{i=1}^{I_{\max }}}{m^{\prime }}_i^G{n^{\prime }}_i^G{V}_i^G}{{\displaystyle \sum _{i=1}^{I_{\max }}}{m^{\prime }}_i^G{n^{\prime }}_i^G}$$(25)

Consequently, the mean ganglion size pertaining to the population of stranded ganglia, ${V^{\prime }}_m^{SG}$, is estimated as $${V^{\prime }}_m^{SG}=\frac{{\displaystyle \sum _{i=1}^{I_{\max }}}\left(1-{m^{\prime }}_i^G\right){n^{\prime }}_i^G{V}_i^G}{{\displaystyle \sum _{i=1}^{I_{\max }}}\left(1-{m^{\prime }}_i^G\right){m^{\prime }}_i^G{n^{\prime }}_i^G}$$(26)

3.2 Interstitial layout of the flow in terms of oil disconnectedness and o/w surface areas

The o/w surface areas are computed per unit volume of porous medium (p.u.v.p.m.), ${{{\displaystyle \tilde{A}}}^{\prime }}_{ow}$, to provide a rough indication of the structure of the entire flow, or per unit volume (p.u.v.) of the disconnected − oil flow domains, ${{{\displaystyle \tilde{A}}}^{\prime }}_{ow,DOF}$, e.g. p.u.v. of GCD, ${{{\displaystyle \tilde{A}}}^{\prime }}_{ow,G}$, or p.u.v. of DTF domain, ${{{\displaystyle \tilde{A}}}^{\prime }}_{ow,D}$. Still more information can be retrieved by accounting the specific (p.u.v.) o/w surface area pertaining to mobilized disconnected-oil blobs, ${{{\displaystyle \tilde{A}}}^{\prime }}_{ow}^M$, these being ganglia, ${{{\displaystyle \tilde{A}}}^{\prime }}_{ow}^{MG}$, and/or droplets, ${{{\displaystyle \tilde{A}}}^{\prime }}_{ow}^{MD}$, so that ${{{\displaystyle \tilde{A}}}^{\prime }}_{ow}={{{\displaystyle \tilde{A}}}^{\prime }}_{ow}^M+{{{\displaystyle \tilde{A}}}^{\prime }}_{ow}^S$ and ${{{\displaystyle \tilde{A}}}^{\prime }}_{ow}^M={{{\displaystyle \tilde{A}}}^{\prime }}_{ow}^{MG}+{{{\displaystyle \tilde{A}}}^{\prime }}_{ow}^{MD}$, or pertaining to stranded oil blobs, ${{{\displaystyle \tilde{A}}}^{\prime }}_{ow}^s={{{\displaystyle \tilde{A}}}^{\prime }}_{ow}^{SG}+{{{\displaystyle \tilde{A}}}^{\prime }}_{ow}^{SD}$. 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 ${{\displaystyle \tilde{A}}}_{ow,i}^{mMG}$. 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_{ow}^{rec}$ and $A_{ow}^{adv}$, 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, $${{\displaystyle \tilde{A}}}_{ow}^{mMG}={{\displaystyle \tilde{A}}}_{ow}^{rec}+{{\displaystyle \tilde{A}}}_{ow}^{adv}$$(27)

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, $${{\displaystyle \tilde{A}}}_{ow,i}^{sMG}=2i\left(N_d-1\right){{\displaystyle \tilde{A}}}_{ow}^{throat}$$(28)

Total mean surface are of the o/w interface of an i-class mobilized ganglion $${{\displaystyle \tilde{A}}}_{ow,i}^{MG}={{\displaystyle \tilde{A}}}_{ow}^{mMG}+{{\displaystyle \tilde{A}}}_{ow,i}^{sMG}$$(29)

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, $${{\displaystyle \tilde{A}}}_{ow,i}^{SG}={{\displaystyle \tilde{A}}}_{ow}^{chamber}+{{\displaystyle \tilde{A}}}_{ow}^{throat}+2i\left(N_d-1\right){{\displaystyle \tilde{A}}}_{ow}^{throat}$$(30)

In the previous expressions, ${{\displaystyle \tilde{A}}}_{ow}^{chamber},{{\displaystyle \tilde{A}}}_{ow}^{throat}$, 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^{\prime }}_{owMMG},{A^{\prime }}_{owSMG},{A^{\prime }}_{owMG}$, ${A^{\prime }}_{owSG},{A^{\prime }}_{owG}$ as well as for o/w interfaces in the DTF domain, A_owMD, A_owSD, A_owD. In all relevant expressions, the index notation consistently determines the physical quantities as follows: G – ganglia, D – Droplets, (ow)MMG/SMG – mobilized/stranded interface of mobilized ganglia, (ow)MD/SD – mobilized/stranded interface of droplets, (ow)G – total interface area of ganglia etc. In particular,  the mobilized o/w surface area of mobilized ganglia p.u.v. of GCD is estimated by $$\frac{1}{{\displaystyle \tilde{\ell }}}{A^{\prime }}_{owMMG}={\displaystyle \sum _{i=1}^{I_{\max }}}{m^{\prime }}_i^G{{{\displaystyle \tilde{n}}}^{\prime }}_i^G{{\displaystyle \tilde{A}}}_{ow,i}^{mMG}={\displaystyle \sum _{i=1}^{I_{\max }}}{m^{\prime }}_i^G{{{\displaystyle \tilde{n}}}^{\prime }}_i^G\left({{\displaystyle \tilde{A}}}_{ow}^{rec}+{{\displaystyle \tilde{A}}}_{ow}^{adv}\right)=\frac{1}{{\displaystyle \tilde{\ell }}}\frac{1}{2N_d\epsilon }\left(A_{ow}^{rec}+A_{ow}^{adv}\right){\displaystyle \sum _{i=1}^{I_{\max }}}{m^{\prime }}_i^G{n^{\prime }}_i^G$$(31) the stranded o/w surface area of mobilized ganglia p.u.v. of GCD is $$\frac{1}{{\displaystyle \tilde{\ell }}}{A^{\prime }}_{owSMG}={\displaystyle \sum _{i=1}^{I_{\max }}}{m^{\prime }}_i^G{{{\displaystyle \tilde{n}}}^{\prime }}_i^G{{\displaystyle \tilde{A}}}_{ow,i}^{sMG}=\frac{1}{{\displaystyle \tilde{\ell }}}\frac{\left(N_d-1\right)}{N_d\epsilon }{A}_{ow}^{throat}{\displaystyle \sum _{i=1}^{I_{\max }}}i{m^{\prime }}_i^G{n^{\prime }}_i^G$$(32)the total o/w surface area of mobilized ganglia p.u.v. of GCD is $$\frac{1}{{\displaystyle \tilde{\ell }}}{A^{\prime }}_{owMG}={\displaystyle \sum _{i=1}^{I_{\max }}}{m^{\prime }}_i^G{{\displaystyle \widetilde{n^{\prime }}}}_i^G{{\displaystyle \tilde{A}}}_{ow,i}^{MG}=\frac{1}{{\displaystyle \tilde{\ell }}}\frac{1}{2N_d\epsilon }{\displaystyle \sum _{i=1}^{I_{\max }}}{m^{\prime }}_i^G{n^{\prime }}_i^G{A}_{ow,i}^{MG}$$(33)

Similarly, the o/w surface area of stranded ganglia p.u.v. of GCD, is $$\frac{1}{{\displaystyle \tilde{\ell }}}{A^{\prime }}_{owSG}={\displaystyle \sum _{i=1}^{I_{\max }}}\left(1-{m^{\prime }}_i^G\right){{{\displaystyle \tilde{n}}}^{\prime }}_i^G{{\displaystyle \tilde{A}}}_{ow,i}^{SG}=\frac{1}{{\displaystyle \tilde{\ell }}}\frac{1}{2N_d\epsilon }{\displaystyle \sum _{i=1}^{I_{\max }}}\left(1-{m^{\prime }}_i^G\right){n^{\prime }}_i^G{A}_{ow,i}^{SG}$$(34)

The stranded o/w surface area of ganglia (mobilized and stranded) p.u.v. of GCD, is $$\frac{1}{{\displaystyle \tilde{\ell }}}\left({A^{\prime }}_{owSMG}+{A^{\prime }}_{owSG}\right)=\frac{1}{{\displaystyle \tilde{\ell }}}\frac{1}{N_d\epsilon }\left[\left(N_d-1\right)A_{ow}^{throat}{\displaystyle \sum _{i=1}^{I_{\max }}}i{m^{\prime }}_i^G{n^{\prime }}_i^G+\frac{1}{2}{\displaystyle \sum _{i=1}^{I_{\max }}}\left(1-{m^{\prime }}_i^G\right){n^{\prime }}_i^G{A}_{ow,i}^{SG}\right]$$(35)

The total o/w surface area of ganglia (mobilized and stranded) p.u.v. of GCD, is $$\frac{1}{{\displaystyle \tilde{\ell }}}{A^{\prime }}_{owG}=\frac{1}{{\displaystyle \tilde{\ell }}}\left({A^{\prime }}_{owMMG}+{A^{\prime }}_{owSMG}+{A^{\prime }}_{owSG}\right)$$(36)

Specific (p.u.v.p.m.) surface areas of o/w-interface in the DTF region

The average (mean) oil drop volume is $$V_{{}^D}{{\displaystyle \tilde{\ell }}}^3={{\displaystyle \tilde{\ell }}}^3{\displaystyle \sum _{j=1}^5}{f}_j{V}_{D,j}={{\displaystyle \tilde{\ell }}}^3{\displaystyle \sum _{j=1}^5}{f}_j\pi {\left(a{D}_{Tj}\right)}^3/6$$(37)

The average number of oil drops p.u.v. of DTF domain is $${{{\displaystyle \tilde{N}}}^{\prime }}_{{}^D}=\frac{1}{{{\displaystyle \tilde{\ell }}}^3}{S^{\prime }}_{{}^{o,D}}\frac{\epsilon }{V_{{}^D}}$$(38)

The o/w surface area of mobilized drops p.u.v. of DTF domain is $$\frac{1}{{\displaystyle \tilde{\ell }}}{A^{\prime }}_{owMD}={m^{\prime }}_{{}^D}{{{\displaystyle \tilde{N}}}^{\prime }}_{{}^D}{A}_{ow}^{mD}$$(39) while the o/w surface area of stranded drops p.u.v. of DTF domain is $$\frac{1}{{\displaystyle \tilde{\ell }}}{A^{\prime }}_{owSD}=\left(1-{m^{\prime }}_{{}^D}\right){{{\displaystyle \tilde{N}}}^{\prime }}_{{}^D}{A}_{ow}^{sD}$$(40)where $A_{ow}^{mD}$ and $A_{ow}^{sD}$ 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_{ow}^{mD}$, 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_{ow}^{sD}$, 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 $$\frac{1}{{\displaystyle \tilde{\ell }}}{A^{\prime }}_{owD}=\frac{1}{{\displaystyle \tilde{\ell }}}{A^{\prime }}_{owMD}+\frac{1}{{\displaystyle \tilde{\ell }}}{A^{\prime }}_{owSD}$$(41)

The total o/w surface area p.u.v. of DOF domain is $$\frac{1}{{\displaystyle \tilde{\ell }}}{A^{\prime }}_{owDOF}=\frac{1}{{\displaystyle \tilde{\ell }}}{\omega }^{\prime }{A^{\prime }}_{owG}+\frac{1}{{\displaystyle \tilde{\ell }}}\left(1-{\omega }^{\prime }\right){A^{\prime }}_{owD}$$(42)

The total o/w surface area p.u.v. of the whole pore network is $$\frac{1}{{\displaystyle \tilde{\ell }}}{A^{\prime }}_{ow}=\frac{1}{{\displaystyle \tilde{\ell }}}\left(1-{\beta }^{\prime }\right){A^{\prime }}_{owDOF}=\frac{1}{{\displaystyle \tilde{\ell }}}\left(1-{\beta }^{\prime }\right)\left[{\omega }^{\prime }{A^{\prime }}_{owG}+\left(1-{\omega }^{\prime }\right){A^{\prime }}_{owD}\right]$$(43)

Expected, ensemble average values of the reduced, specific o/w surface area in the DOF domains, for different values of the capillary number, the flow rate ratio and the viscosity ratio, are presented in Section 4.

Coefficient (fraction) of oil fragmentation

A measure of the extent of oil disconnection is provided by the coefficient of oil fragmentation, ${f^{\prime }}_{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 thoat size classes. The minimum size of an oil-droplet is, $$V_{D,1}{{\displaystyle \tilde{\ell }}}^3={{\displaystyle \tilde{\ell }}}^3\pi {\left(D_{T1}\right)}^3/6$$(44) whereas the corresponding o/w surface area is $${{\displaystyle \tilde{\ell }}}^2\pi {\left(D_{T1}\right)}^2$$(45)

Now, consider the number density of the oil drops that could be formed if all of the oil in place had shuttered into oil droplets of minimum size, ${\tilde{\mathit{V}}}_{\mathit{D},1}{{\displaystyle \ell }}^3$, $${N^{\prime }}_{crit}^D=\frac{{{\displaystyle \tilde{\ell }}}^3\epsilon \left(1-{S^{\prime }}_w\right)}{{{\displaystyle \tilde{\ell }}}^3{V}_{D,1}}=\left(1-{S^{\prime }}_w\right)\frac{\epsilon }{V_{D,1}}$$(46)

The respective specific o/w surface area, ${{{\displaystyle \tilde{A}}}^{\prime }}_{ow,crit}$, would comprise spherical surfaces pertaining to the mimimum size drops, of diameter equal to the diameter of the minimum throat, ${\tilde{\mathrm{D}}}_{h1}$, therefore, $$\frac{1}{{{\displaystyle \tilde{\ell }}}^3}{{{\displaystyle \tilde{A}}}^{\prime }}_{ow,crit}=\frac{1}{{{\displaystyle \tilde{\ell }}}^3}{N^{\prime }}_{crit}^D{{\displaystyle \tilde{\ell }}}^2\pi {\left(D_{h1}\right)}^2=\frac{1}{{\displaystyle \tilde{\ell }}}\left(1-{S^{\prime }}_w\right)\frac{\epsilon }{V_{D,1}}\pi {\left(D_{h1}\right)}^2=\frac{1}{{\displaystyle \tilde{\ell }}}{A^{\prime }}_{ow,crit}$$(47) The coefficient (or, fraction) of oil fragmentation is defined as $${f^{\prime }}_{OF}=\frac{{A^{\prime }}_{ow}}{{A^{\prime }}_{ow,crit}}$$(48)

It is straightforward to derive its limiting values as $0\le {f^{\prime }}_{OF}\le 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^{\prime }}_{OF}$ are presented in Section 4.

3.3 Interstitial flow structure in terms of transport properties

To map the flow structure in terms of transport properties and not just the flow layout, one should account for dynamic variables that can efficiently describe the transport characteristics of the flow. In particular, ganglion and droplet velocity distribution, fractional flows of connected (CPF) and 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, $${{\eta }^{\prime }}_{o,C}=\frac{{{{\displaystyle \tilde{q}}}^{\prime }}_{o,C}}{{{\displaystyle \tilde{q}}}_o}={\beta }^{\prime }{U^{\prime }}_{o,C}$$ $${{\eta }^{\prime }}_{o,G}=\frac{{{{\displaystyle \tilde{q}}}^{\prime }}_{o,G}}{{{\displaystyle \tilde{q}}}_o}=\left(1-{\beta }^{\prime }\right){\omega }^{\prime }{U^{\prime }}_{o,G}$$ $${{\eta }^{\prime }}_{o,D}=\frac{{{{\displaystyle \tilde{q}}}^{\prime }}_{o,D}}{{{\displaystyle \tilde{q}}}_o}=\left(1-{\beta }^{\prime }\right)\left(1-{\omega }^{\prime }\right){U^{\prime }}_{o,D}$$(49) where ${U^{\prime }}_{o,G}$ and ${U^{\prime }}_{o,D}$ 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^{\prime }}_{o,DOF}$, as the former comprise an inseparable mixture within DOF. They are given by $${U^{\prime }}_{o,G}={U^{\prime }}_{o,D}={U^{\prime }}_{o,DOF}=\frac{1-{\beta }^{\prime }{U^{\prime }}_{o,C}}{1-{\beta }^{\prime }}$$(50)where ${U^{\prime }}_{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, η_o,C, GD, η_o,G, and DTF, η_o,D, for different values of the capillary number, the flow rate ratio and the viscosity ratio, are presented in Section 4.

Superficial velocity of o/w surface in the GD domain

The superficial velocity of the o/w interface of the entire swarm of ganglia,${{{\displaystyle \tilde{U}}}^{\prime }}_{ow,G}$, may be estimated as the total momentum of ganglia o/w surface divided by the total o/w surface area, $${{{\displaystyle \tilde{U}}}^{\prime }}_{ow,G}=\frac{{\displaystyle \sum _{i=1}^{I_{\max }}}{m^{\prime }}_i^G{{{\displaystyle \tilde{n}}}^{\prime }}_i^G{{\displaystyle \tilde{A}}}_{ow,i}^{mMG}{{{\displaystyle \tilde{u}}}^{\prime }}_i^G}{{\displaystyle \sum _{i=1}^{I_{\max }}}{{{\displaystyle \tilde{n}}}^{\prime }}_i^G{{\displaystyle \tilde{A}}}_{ow,i}^G}=\frac{{\displaystyle \sum _{i=1}^{I_{\max }}}{m^{\prime }}_i^G{{{\displaystyle \tilde{n}}}^{\prime }}_i^G{{\displaystyle \tilde{A}}}_{ow,i}^{mMG}\left({u^{\prime }}_i^G{{\displaystyle \tilde{U}}}_{{}^w}\right)}{{\displaystyle \sum _{i=1}^{I_{\max }}}{{{\displaystyle \tilde{n}}}^{\prime }}_i^G{{\displaystyle \tilde{A}}}_{ow,i}^G}={{\displaystyle \tilde{U}}}_{{}^w}\frac{{\displaystyle \sum _{i=1}^{I_{\max }}}{m^{\prime }}_i^G{n^{\prime }}_i^G{A}_{ow,i}^{mMG}{u^{\prime }}_i^G}{{A^{\prime }}_{ow,G}}$$(51) and then, reduced by the macroscopic velocity of oil, it becomes: $${U^{\prime }}_{ow,G}=\frac{1}{{{\displaystyle \tilde{U}}}_o}{{{\displaystyle \tilde{U}}}^{\prime }}_{ow,G}=\frac{1}{r}\frac{{\displaystyle \sum _{i=1}^{I_{\max }}}{m^{\prime }}_i^G{n^{\prime }}_i^G{A}_{ow,i}^{mMG}{u^{\prime }}_i^G}{{A^{\prime }}_{ow,G}}$$(52)

Superficial velocity of o/w surface in the DTF domain

The reduced superficial velocity of the o/w surface in the DTF domain is readily estimated as $${U^{\prime }}_{ow,mD}=A_{ow}^{mD}{U^{\prime }}_{o,D}$$(53)

It is worth noting the difference between the expressions in Equations (52) and (53). The summation in Equation (52) is to account the different momenta of the ganglion size classes migrating downstream at different velocities (to each ganglion size class corresponds a size class velocity), whereas in the indistinguishable DTF unit cells all of the drops that are mobilized migrate downstream at the uniform superficial velocity, ${U^{\prime }}_{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 $${U^{\prime }}_{ow,DOF}={\omega }^{\prime }{U^{\prime }}_{ow,G}+\left(1-{\omega }^{\prime }\right){U^{\prime }}_{ow,D}$$(54)

Fractional o/w surface transport through GD and DTF

The fractional transport of the o/w interfaces through GD and DTF, ${{\xi }^{\prime }}_{ow,G}$ and ${{\xi }^{\prime }}_{ow,D}$, may be estimated as the momentum of ganglia o/w surface divided by the total o/w surface area, $${{\xi }^{\prime }}_{ow,G}==\frac{{\omega }^{\prime }{U^{\prime }}_{owG}}{{\omega }^{\prime }{U^{\prime }}_{owG}+\left(1-{\omega }^{\prime }\right){U^{\prime }}_{owD}}\text{ and }{{\xi }^{\prime }}_{ow,D}=1-{{\xi }^{\prime }}_{ow,G}$$(55)

Expected, ensemble average values of ${U^{\prime }}_{ow,DOF}$, ${{\xi }^{\prime }}_{ow,G}$ and ${{\xi }^{\prime }}_{ow,D}$ for different values of the capillary number, the flow rate ratio and the viscosity ratio, are presented in Section 4.

4 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, ${{\displaystyle \tilde{\mu }}}_w=0.94\times {10}^{-3}\text{Pa}\text{s}$, oil/water interfacial tension ${{\displaystyle \tilde{\gamma }}}_{ow}=25\times {10}^{-3}\text{N}/\text{m}$, advancing (A) and receding (R) oil-water/p.m. contact angles, ${\theta }_A^0=43\mathrm{deg}$ and ${\theta }_R^0=37\mathrm{deg}$. Five typical systems were examined: favorable, κ = 0.33, 0.67, indifferent κ = 1.00 and unfavorable, κ = 1.50 and 3.00, viscosity ratio values were examined.

The simulations covered a rectangular domain logCa   ×   logr, from (logCa, logr) = (− 8, − 2) to (logCa, logr) = (− 3, + 2). The domain was covered in 9 successive steps per order of capillary number magnitude, i.e. Ca = 1,  ...,  9 × 10^−j,  j = − 8,   ...,   − 4 (9 × 5 =45 divisions in the logCa range), and in 21 successive steps of logr = 0.2 in the logr range. Therefore, 45 × 21 = 945 simulations, pertaining to different flow conditions, were run for each of the 5 different values of the viscosity ratio. Each diagram in Figure 10 and in 12–19 contain the points that depict the average values of physical quantities (dependent variables) describing the interstitial structure of the flow for the corresponding set of flow conditions (Ca, r). Each diagram in these figures should normally contain 945 points. Nevertheless, this number is not reached as there are many flow conditions for which no solutions could be detected, i.e. 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 κ = 0.33, κ = 1.00 and κ = 3.00 are presented. The flow structure maps of the other two systems (κ = 0.67 and 1.50) fall in between.

The afore mentioned diagrams (Fig. 10, 12–19) share common characteristics. The abscissa axis measures the logarithm of the flow rate ratio. In every diagram, iso-capillary 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⁻⁸ ≤ Ca ≤ 10⁻⁴, while a dashed line arrow is for changes over 4 × 10⁻⁷ ≤ Ca ≤ 10⁻⁴ (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⁻⁴

• For equal viscosity systems (κ = 1), the DeProF algorithm failed to detect any physically admissible solutions in the low Ca regime; specifically no solution could be detected for Ca < 4 × 10⁻⁷. Presumably, viscosity disparity allows the process more degrees of freedom in adjusting itself within a richer ensemble of equiprobable physically admissible flow configurations (see discussion, 4.10).

The following Sections 4.1–4.8 provide an overall (basic) description of the furnished diagrams. Extended discussion of the particular trends in the flow arrangement/structure, as observed across flow conditions and system configurations, are provided in Section 4.9.

4.1 Reduced pressure gradient

Expected, ensemble average values of the reduced pressure gradient, x, Equation (5), for different values of the capillary number, Ca, and the flow rate ratio, r, are presented in the diagrams of Figure 10 for different viscosity ratio systems, κ.

Observations of single phase flows within pore networks confirm that the macroscale pressure gradient scales linearly with the superficial fluid velocity governed by Darcy's law. This seems to be a quite trustworthy modeling consideration also in the case of 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 attributed to differences associated to (a) dimensionality of the pertinent variables (b) measurements pertaining to different flow conditions, (c) dimensionality of the o/w/pm system properties. The DeProF model predictions bridge these discrepancies and show the potential for deriving a universal power law between appropriate dimensionless variables.

4.2 Typical, physically admissible flow configurations

Indicative physically admissible flow configurations described by the values of the FAV, $({S^{\prime }}_w,{\beta }^{\prime },{\omega }^{\prime })$, and the associated values of the ganglion size distribution parameters, $({n^{\prime }}_1,{n^{\prime }}_2,{\zeta }^{\prime })$, for different flow conditions, are presented in the diagrams of Figure 11. ‘Moderate flow arrangements’ are characterized by physically admissible flow arrangement values $({S^{\prime }}_w,{\beta }^{\prime },{\omega }^{\prime })$ close to the ensemble average, (S_w, β, ω) and are presented in diagrams (a1)–(c1). ‘Extreme flow arrangements’ are characterized by physically admissible flow arrangement values $({S^{\prime }}_w,{\beta }^{\prime },{\omega }^{\prime })$ 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, β, the volume fraction of GD within the DOF, ω, for different values of the capillary number, Ca, and the flow rate ratio, r, for different viscosity ratio systems, κ, are presented in the diagrams of Figure 12. As the viscosity ratio is increased, the onset and build-up of connected oil pathways (β) are observable from lower Ca values.

4.4 Reduced mechanical power dissipation

The reduced rate of mechanical energy dissipation within the process, W, is defined as (Valavanides et al., 2016) $$W=\frac{{\displaystyle \tilde{W}}}{{{\displaystyle \tilde{W}}}^{1\mathrm{\Phi }}}={\displaystyle \tilde{W}}\frac{{\displaystyle \tilde{k}}{{\displaystyle \tilde{\mu }}}_w}{{\left({{\displaystyle \tilde{\gamma }}}_{ow}Ca\right)}^2}=\kappa r\frac{1}{k_{ro}}(r+1)=\frac{1}{k_{rw}}(r+1)=\left[\beta r{U}_{o,C}+(1-\beta )\left(r{U}_{o,G}+U_{w,G}\right)\right]x$$(56) where ${{\displaystyle \tilde{W}}}^{1\mathrm{\Phi }}={\left({{\displaystyle \tilde{\gamma }}}_{ow}Ca\right)}^2/\left({\displaystyle \tilde{k}}{{\displaystyle \tilde{\mu }}}_w\right)$ is the mechanical power dissipation for the equivalent one-phase flow of water with superficial velocity, ${{\displaystyle \tilde{U}}}_w$. The mechanical power dissipated in Connected Pathway Flow (CPF), Ganglion
Dynamics (GD) and Drop-Traffic Flow (DTF) are given respectively by $$W_C=\beta r{U}_{o,C}x$$ $$W_G=\omega (1-\beta )\left(r{U}_{o,G}+U_{w,G}\right)x$$(57) $$W_D=(1-\omega )(1-\beta )\left(r{U}_{o,G}+U_{w,G}\right)x$$

Expected values of the reduced mechanical power dissipation, W, and its partition into the individual interstitial, constituent flows, for different values of the capillary number, Ca, and the flow rate ratio, r, for different viscosity ratio systems, κ_, are presented in the diagrams of Figure 13.

4.5 Mean ganglion size (overall and mobilized)

Expected, values of the mean ganglion sizes, averaging Equations (24) and (25) over the ensemble of physically admissible solutions, for different values of the capillary number, Ca, and viscosity ratio, κ, are presented in the diagrams of Figure 14. The mean size pertaining to the entire ganglia population (mobilized and stranded), $V_m^G$, is presented in (a), whereas the mean size of the mobilized ganglia,$V_m^{mG}$, in (b).

4.6 Interfacial areas and oil disconnectedness

Expected, ensemble average values of the reduced o/w interface surface area for different values of the capillary number, Ca, flow rate ratio, r, and viscosity ratio, $\kappa ={{\displaystyle \tilde{\mu }}}_o/{{\displaystyle \tilde{\mu }}}_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, ${{\displaystyle \tilde{A}}}_{ow}^m$, whether these are ganglia, ${{\displaystyle \tilde{A}}}_{ow}^{mG}$, or droplets, ${{\displaystyle \tilde{A}}}_{ow}^{mD}$, so that ${{\displaystyle \tilde{A}}}_{ow}={{\displaystyle \tilde{A}}}_{ow}^m+{{\displaystyle \tilde{A}}}_{ow}^s$ and ${{\displaystyle \tilde{A}}}_{ow}^m={{\displaystyle \tilde{A}}}_{ow}^{mG}+{{\displaystyle \tilde{A}}}_{ow}^{mD}$. 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.

The extent of oil disconnectedness can be accounted for by the coefficient of oil fragmentation, f_OF, Equation (48). Expected values of f_OF, for different values of the capillary number, Ca, flow rate ratio, r, and viscosity ratio, κ, are presented in the diagrams of Figure 16. In particular, f_OF can be used to assess the unattained margin of o/w surface area density. It can be considered as a reduced A_ow, normalized by the maximum value of A_ow that would be attained if all of the oil would shutter in droplets of the smallest possible size. In that context, f_OF could also be called ‘interfacial saturation’. As the flow rate ratio, r, increases, f_OF drops significantly, indicating that a part of DOF transforms (reorganizes) into 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.

4.7 Dynamic variables describing the interstitial transport properties of the process

Ensemble average values of the fractional flow of oil through CPF, η_o,C, GD, η_o,G, and DTF, η_o,D, Equations (49), are presented in the diagrams of Figure 17, for different values of the capillary number, the flow rate ratio and the viscosity ratio.

In all cases examined, the contribution of disconnected-oil 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, ξ_ow,D, and GD, ξ_ow,G, Equations (54) and (55), are presented in the diagrams of Figure 18, for different values of the capillary number, the flow rate ratio and the viscosity ratio.

In all cases, the contribution of transport of o/w interface through DTF, prevails over GD; in particular, across the low-end domain of flow rate ratio, r, the o/w interface transport takes place almost exclusively through DTF. This is attributed to the large volumetric density of o/w interfaces in DTF, Figure 15, in combination to the fractional flow of disconnected oil through DTF, Figure 17, across the entire domain of flow conditions.

4.8 Energy efficiency of 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, Eqs. (6)–(8)] as $$f_{EU}(Ca,r)=\frac{{{\displaystyle \tilde{U}}}_o}{{{\displaystyle \tilde{U}}}_w}\frac{{\left({{\displaystyle \tilde{\gamma }}}_{ow}Ca\right)}^2}{{\displaystyle \tilde{W}}{\displaystyle \tilde{k}}{{\displaystyle \tilde{\mu }}}_w}=\frac{r}{W(Ca,r)}$$(58)

Simulations implementing the DeProF model, suggest that a set of flow conditions exist for which the sought process attains locally maximum energy efficiency, see Figure 19. Such conditions, define a locus of optimum operating conditions or critical flow conditions, r*(Ca), for each oil-water-p.m. system (Valavanides, 2014 and Valavanides, 2018). The locally maximum value in f_EU,ij (Ca_i, r_j) detected in the simulations per fixed Ca_i, corresponds to a particular flow rate ratio value, $r_i^{*}$.The sets $\left\{C{a}_i,r_i^{*}\right\}$ 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 viscosity-dominated. 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’), $f_{EU\infty }^{*}=1/{\left(1+\sqrt{\kappa }\right)}^2$, while the vertical, dashed line in red indicates the corresponding flow rate ratio, $r_{\infty }^{*}=1/\sqrt{\kappa }$. 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 steady-state two-phase flows in real porous media, remaining in latency until indicated by DeProF model predictions (Valavanides and Payatakes, 2003). Recently, in order to develop a theoretical framework based on statistical thermodynamics that would justify the existence of a locus optimum operating conditions, Valavanides and Daras (2016) proposed a model for accounting the number of the process microstates.

4.9 Discussion of results

The diagrams presented in Figures 10–19 indicate that distinct mutations can be observed. Some mutations are obvious and measurable at the macroscopic scale, others remain in latency within the interstitial flow configurations and, in practice, they are very difficult to measure. Nevertheless, the DeProF model has the capability to predict what could be revealed in a systematic laboratory study, would the necessary equipment and resources be available.

The behavior of the 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⁻⁷ for medium to 10⁻⁴ for very low Ca). Getting deeper in the interstitial flow configuration, a mutation can be revealed as the flow rate ratio is increased from relatively low to higher values. Specifically,

Low Ca − low r − At the low end of the flow rate ratio domain, the prevailing flow pattern is (DOF) as practically no (CPF), β → 0, can be detected in the diagrams of row (b), Figure 12. Delving deeper in DOF, DTF occupies a significantly larger − over GD − fraction of the available p.m. cross-section, as is indicated by the rather low ω and high (1 − ω) values predicted in the diagrams of Figure 12c. Still for relatively low r values, DTF is also the prevailing flow pattern as is indicated in the oil fractional flows diagrams, Figure 17, where CPF is practically 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_m^{mG}\approx 10$ [diagrams in row (b)], large enough, if compared to the mean size of all ganglia,$V_m^G\approx 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 − β) → 1, and borderline CPF β → 0 for −2 < logr < − 1, to mixed DOF and CPF for −1 < logr < 0, see Figure 12b. Delving deeper in DOF, DTF occupies a predominant fraction of the available p.m. cross-section − over GD, as is indicated by the low ω and high (1 − ω) values on the diagrams of Figure 12c across the whole range of flow rate ratio values (−2 < logr < 0). Albeit CPF does not occupy a significant volume fraction (topping up to 0.22 for κ = 0.33), its contribution in the oil transport is dramatic as it suddenly 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 κ = 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 $\left(V_m^G\approx 1.0\right)$, see Figure 14a. The mean size of mobilized ganglia is also small $\left(2<V_m^{mG}<3\right)$, Figure 14b, because of the large pressure gradients. Nevertheless, the low volumetric occupation of GD in combination with the low migration velocity of small ganglia (as the driving force − the externally induced net pressure drop − is small) result in small fractional flow through GD, as attested in Figure 17b. The fractional flow of oil through DTF (the GD complementary constituent in DOF) is significant; nevertheless, it drops dramatically as CPF takes over, as attested in the diagrams in Figure 17c and a respectively. These statements conform with the predicted values of the o/w surface area densities, Figure 15, of the coefficient of oil fragmentation, Figure 16, and the obvious prevalence of DTF as the main transporter of o/w interface over GD, as shown in Figure 18, diagrams in row (b) against those in row (c).

Volume fractions against fractional flows

The previous analysis leads to the issue of whether using volume fractions in describing the flow, is more specific than using fractional flows. If we focus on the description of connected oil flow (the CPF domain), we may observe significant discrepancies between volume fractions and fractional flows − already identified in the previous paragraphs. In general, for the same flow conditions and viscosity ratios, the predicted values of the fractional flow rates, η_o,C, are larger than the corresponding values of the volume fractions, β. This is expected if we consider that the absence of o/w interfaces in the CPF domain does not impede the flow (as the the bulk viscosity does). Let us compare the diagrams of the volume fraction of connected-oil pathways, β, Figure 12b, with those of the fractional flow of oil through CPF, η_o,C, Figure 17b. Indicative values are tabulated for typical flow conditions and viscosity ratios, see Table 2. The discrepancies are pronounced in the low Ca / low r regimes and they tend to 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.

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 iso-capillary 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 (κ = 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).

Systematic mutations have been revealed within the flow structure 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. Τhis 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 (κ ≠ 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 non-wetting 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.

Acknowledgment

This research work has been co-funded by the European Union (European Social Fund) and Greek national resources under the framework of the “Archimedes III: Funding of Research Groups in TEI of Athens” project of the “Education and Lifelong Learning” Operational Program (MlS379389).

Dr. Olga Vizika, IFP Energies Nouvelles, is greatly acknowledged for discussions on the subject. Mrs. Evangelia Ponta, TEI Athens, is greatly acknowledged for assisting with the simulations. Valuable remarks and constructive comments from anonymous referees are also gratefully acknowledged.

References

All Tables

All Figures