Flow Rate Dependency of Steady-State Two-Phase Flows in Pore Networks: Universal, Relative Permeability Scaling Function and System-Characteristic Invariants

The phenomenology of steady-state two-phase flow in porous media is conventionally recorded by the relative permeability diagrams in terms of saturation. Yet, theoretical, numerical and laboratory studies of flow in artificial pore network models and natural porous media have revealed a significant dependency on the flow rates—especially when the flow regime is capillary to capillary/viscous and part of the disconnected non-wetting phase remains mobile. These studies suggest that relative permeability models should incorporate the functional dependence on flow intensities. In the present work, a systematic dependence of the pressure gradient and the relative permeabilities on flow rate intensity is revealed. It is based on extensive simulations of steady-state, fully developed, two-phase flows within a typical 3D model pore network, implementing the DeProF mechanistic–stochastic model algorithm. Simulations were performed across flow conditions spanning 5 orders of magnitude, both in the capillary number, Ca, and the flow rate ratio, r, and for different favorable /unfavorable viscosity ratio fluid systems. The systematic, flow rate dependency of the relative permeabilities can be described analytically by a universal scaling function along the entire domain of the independent variables of the process, Ca and r. This universal scaling comprises a kernel function of the capillary number, Ca, that describes the asymmetric effects of capillarity across the entire flow regime—from capillarity-dominated to mixed capillarity/viscosity- to viscosity-dominated flows. It is shown that the kernel function, as well as the locus of the cross-over relative permeability values, are single-variable functions of the capillary number; they are both identified as viscosity ratio invariants of the system. Both invariants can be correlated with the structure of the pore network, through a function of Ca. Consequently, the correlation is associated with the wettability characteristics of the system. Among the potential applications, the proposed, universal, flow rate dependency scaling laws are the improvement of core analysis and dynamic rock-typing protocols, as well as integration into field-scale simulators or associated machine learning interventions for improved specificity/accuracy.


Introduction
Τwo-phase flow in porous media is a physical process whereby two phases flow simultaneously within a porous medium.A porous medium (PM) comprises a complex network of interconnected pores of sizes in the range of microns-to-millimeters.When the flow is immiscible, i.e., the two phases do not mix, one of the phases is wetting the interstitial surface of the pore network against the other, non-wetting phase.
Applications of two-phase flow in porous media are numerous within the industry, energy and environment sectors, e.g., enhanced oil recovery, CO 2 sequestration and geothermal harvesting, groundwater and soil contamination and subsurface remediation, the operation of multiphase trickle-bed reactors, the operation of proton exchange membrane fuel cells.
The size and time scales of the process span across many orders of magnitude.At the pore scale, the characteristic size is such that the flow takes place at the continuum scale (Stokes flow).Critical flow mechanisms occur within size domains spanning micro-to-millimeter scales.An intrinsic characteristic of immiscible two-phase flow in porous media is the disconnection of the non-wetting phase.Fluid immiscibility, wettability and the complex pore network structure (tortuosity, genus, dimensionality, micro-roughness, surface chemistry, etc.) result in the disconnection of the non-wetting phase into fluidic elements, namely ganglia and droplets, of correspondingly larger or smaller size when compared to the average pore size.Pore scale hydraulic instabilities-instigated by Haines jumps and the corresponding, propagating pressure pulses-trigger detachment of the non-wetting phase, break-up and pinch-off (Payatakes 1982;Berg et al. 2013) and lead to intermittent flows of the fluidic elements of the disconnected non-wetting phase (Reynolds et al. 2017;Spurin et al. 2019) as well as pressure waves and flow instabilities (Hönig et al. 2013;Rücker et al. 2021;Karadimitriou et al. 2023).These mechanisms have their share on increasing the incessant power dissipation in the incoherent motion of fluidic elements.The Young-Laplace forces acting across the N/W menisci separating the fluidic elements of the non-wetting phase from the surrounding wetting phase within the network capillaries are comparable to the Stokes flow, viscosity-induced resistances.As a result, an average pressure gradient builds-up along the flow direction.Averaging is considered over time and space in terms of appropriately sized, representative time periods and volume elements.The contribution of these effects on the macroscopic pressure gradient(s) depends on the interstitial structure of the flow and the interaction, through momentum exchange, between the two fluids and the solid matrix.Viscosity-and capillarity-induced forces have key roles in the momentum exchange.
Depending on the flow conditions, different flow regimes establish.These have been identified both in model pore networks (Avraam and Payatakes 1995;Tsakiroglou et al. 2007) and in natural porous media (Armstrong et al. 2017).In those regimes, a part of the disconnected elements of the non-wetting phase in the form of ganglia and droplets (Payatakes 1982;Tallakstad et al. 2009;Georgiadis et al. 2013;Yiotis et al. 2013;Armstrong et al. 2016) eventually organize to form larger structures at the high-end of the size spectrum, namely connected pathways of the non-wetting phase (Rucker et al. 2015).Similarly, at the low-end of the size spectrum, populations of very small droplets of the nonwetting phase (Datta et al. 2014;Pak et al. 2015) may become so small and form emulsions with the wetting phase (Guillen et al. 2012;Unsal et al. 2016;Rucker et al. 2017;Arab et al 2018).Moving-up from the scale of a few thousands pores into the realm of cm and m (pore-to-core-to-formation), porous medium heterogeneities and discontinuities appear, e.g., cracks, fractures, layering, etc.The collective behavior of the interaction between the two phases and the structure of the pore network has adverse effects on the meter-to-thehundreds-of-meters-size domain.Two-phase flow in natural porous media then becomes a complex system and relevant issues on complexity, emergence, self-organization, etc., need to be considered (Cushman 1997;Faybishenko et al. 2015).So, long as the process may be observed and described at different size scales, it is essential to describe the flow at the representative elementary volume (REV) over the scale associated with the problem of interest.The REV is associated with averaging or up-scaling all interacting flow phenomena occurring within any statistically adequate volume and time frame and, in statistical mechanics terms, it must represent all possible fluidic microstates (Valavanides and Daras 2016;Bedeaux and Kjelstrup 2022).
The industry standard for the analysis and modeling of immiscible, two-phase flows in porous media is built around the concept of relative permeabilities.These have been contrived to extend the classic Darcy's law for one-phase (saturated) flow in porous media, into its fractional form version. Relative permeabilities are calculated from laboratory measurements, either with steady-or unsteady-state methods.In steady-state methods, the two phases are simultaneously injected at a fixed ratio into a core.When the system reaches steady-state conditions, the differential pressures are measured and the relative permeabilities are calculated.The resulting description is phenomenologic.In practice, e.g., in Special Core Analysis (SCAL), an appropriate number of specially selected cores at the size of centi-to-deci-meters are collected from the porous medium formation and then used as a model-or reference-medium to record the effects of simultaneous flow of model fluids (American Petroleum Institute 1998;McPhee et al. 2015).The core is placed tightly in a core holder and fluids are co-injected through the core at controlled flow rates.Measured values of flow rates versus pressure drops are then used to calculate the corresponding values of the relative permeabilities, so as to satisfy the Darcy fractional flow equations [see Eqs.(1), Sect.2.1].This is repeated for different flow conditions (e.g., by imposing different flow rates) in order to derive look-up tables or diagrams that describe the behavior of the examined flow system under immiscible twο-phase flow conditions.In recent research efforts, sophisticated micro-tomographic technologies (Georgiadis et al. 2013;Wildenschild and Sheppard 2013;Fusseis et al. 2014) have been implemented to scan the pore network structure, reveal the saturation profile and the interstitial flow structure along the core (Youssef et al. 2014(Youssef et al. , 2017)), or even to capture micro-events occurring at high resolution domains in terms of time and size scales, e.g., Haines jumps (Berg et al. 2013) or in-situ dynamic alterations of contact angles (Andrew et al. 2015;Singh et al. 2016).Obviously, such improved SCAL techniques require more elaborate interventions and support equipment and are limited to relatively smaller core sizes.The output of SCAL measurements is a core-scale REV model (phenomenologic); the behavior of the core REV model is then incorporated into field-scale flow simulators.
Research efforts have dealt with the problem of correlating relative permeability values to values of saturation, assuming a priori that relative permeabilities are saturationdependent.The assumption of saturation dependency has been established from the early years of oil extraction, and still is the standard approach in the industry until today.The objective is to deliver functional forms for fitting laboratory measured saturation and relative permeability data charts, with scope to incorporate these, together with saturation charts of capillary pressure, into field-scale simulators.Established approaches on fitting relative permeability vs saturation data sets are based on interpolating relative permeability against saturation values, using power law relationships (Corey 1954;Brooks and Corey 1964;Chierici 1984;Honarpour et al. 1986;Lomeland et al. 2005;Lomeland and Ebeltoft 2008;Lomeland 2018).These models attempt to provide appropriate functional expressions for fitting relative permeability values against saturation over the broadest domain possible in flow conditions.The Brooks-Corey and LET (Lomeland-Ebeltoft-Thomas) approaches are most commonly used in fitting data collected during core analysis.Their performance in adequately describing the phenomenology revealed by core analysis, as well as issues associated with sensitivity and uncertainty analysis of relative permeability parameterization, is examined and discussed in a recent publication by Berg et al. (2021).Albeit saturation dependency of relative permeabilities is established in industry, there is skepticism on whether it is the most efficient approach in describing the process.As saturation describes a volumetric state of the system but not a transport state, it brings no definite input to the momentum balance; therefore, it is questionable if it can provide any information on the kinetics of the macroscopic flow.In that context, it is questionable if saturation can be used as an independent variable per se (Valavanides 2018b).The long standing practice in defining the process state through combined relative permeability-capillary pressure-saturation charts does not improve things.The conventional capillary pressure-saturation data sets are obtained under static conditions; even so, they are not unique, as-due to hysteresis-they depend on the direction of the displacement, e.g., drainage or imbibition (Schluter et al. 2016;McClure et al. 2018).
The flow rate dependency of relative permeabilities has been identified as early as the first attempts to reveal the phenomenology of two-phase flow in porous media by Wyckoff and Botset (1936).Nevertheless, a systematic laboratory study across different flow regimes was never followed upon, possibly because of cost and time issues in deploying such an extensive, multi-parametric study and because, at the early days of the oil era, conditions were in consistency-to a certain degree-with saturation dependency.The last 2 decades, the issue of flow dependency is reconsidered, as models based on saturation dependency alone proved to be inefficient in describing the flow at intermediate flow conditions, i.e., when capillarity-induced resistances are comparable to viscous resistances.Theoretical, numerical, and laboratory studies of flow in artificial pore network models and natural porous media (Wildenschild et al. 2001;Nguyen et al. 2006;Aursjo et al. 2014;Sinha et al. 2017;Tsakiroglou et al. 2015;Krause and Benson 2015;Valavanides 2018a;Gao et al. 2020;Yiotis et al. 2019Yiotis et al. , 2021;;Zhang et al. 2021) have identified flow rate dependency to be an underlying, inherent process characteristic, revealing a significant dependency of relative permeabilities on flow rates.
In parallel, research efforts to characterize a porous medium, implementing laboratory measurements under dynamic conditions (flow), namely dynamic rock typing, are currently in development (Mirzaei-Paiaman and Ghanbarian 2021).Studies have also shown the need to establish a comprehensive understanding of the degree and extent of how pertinent system characteristics and process parameters (interfacial tension, wettability, pressure, temperature, etc.) affect the relative permeabilities (Esmaeili et al 2019(Esmaeili et al , 2020;;Kumar et al. 2020;Arab et al. 2018Arab et al. , 2020Arab et al. , 2021).Yet, extracting the flow rate dependency characteristics from physical experiments in the laboratory is not an easy task, so long as the number of essential parameters, e.g., the viscosities of the two fluids, interfacial tension, wettability, do not scale similarly with flow intensity and pore network geometry.That increases significantly the number of independent variables and overloads the protocol of a systematic laboratory study.Moreover, because of end-effects in SCAL methods, there are issues associated with biasing of measurements.Virtual experiments are an alternative approach in dealing with the problem of scaling-up transport properties within a porous medium.Dynamic pore network simulators are an option; nevertheless, the computational time increases exponentially with network size-directly associated with the minimum size of the REV-as well as the degree of disconnection of the non-wetting phase.
Ιn the present work, we attempt to tackle the problem of flow rate dependency.We will see how a generic form of scaling function can be derived in order to describe the dependence of the reduced pressure gradient, x, or, equivalently, the relative permeabilities, on the local flow rate intensities of the wetting and the non-wetting phases.The latter is expressed in reduced form as the capillary number, Ca, and the N/W flow rate ratio, r.
The proposed universal scaling, Eq. ( 18), integrates the contributions of the two main sources of flow resistance, the capillary resistances and the viscous resistances, in decoupled, appropriately reduced format.In particular, the contribution of viscous resistances is accounted by a linear term of the flow rate ratio, r, with coefficient of linear increase the viscosity ratio of the fluid system.Τhe contribution of capillary resistances is accounted by a strongly decreasing function of the capillary number, Ca.This function provides an integrated description of the effects of capillarity related phenomena occurring during the interaction of the interstitial flows within the pore network.It is conceptually identified as the Intrinsic Dynamic Capillary Pressure (IDCP) function of the flow system.The IDCP function can be used as an identity of the fluid system.It is invariant to viscosity ratio.Another characteristic property is the locus of flow conditions whereby equal permeabilities are attained, namely the locus of cross-over relative permeabilities, also invariant to viscosity ratio.Both can be extracted from laboratory or virtual, steady-state co-injection experiments on any given flow system, implementing conventional SCAL techniques.Here, we will use the extensive data output produced from previous, systematic simulations, based on the DeProF model for steady-state, fully developed, two-phase flow in pore networks, Valavanides (2018a).
The novelty of the proposed scaling relays into improvements on the Darcy-scale description of steady-state, fully developed, two-phase flow in pore networks.The current work is complementary to a previous work, implementing an energy efficiency analysis, Valavanides (2018b), for two-phase flow in porous media.In that work, we proved the existence of the locus of critical flow conditions (CFCs), whereby maximum utilization of the input hydraulic power is achieved, as well as the collapsing of all the loci pertaining to different viscosity ratios into a universal locus, pertaining to equal viscosities.Each locus of CFCs is unique; therefore, it can also provide the identity of the flow system, comprising the two fluids and the porous medium, using a different norm (energy efficiency).In the current work, we have implemented fractional flow analysis and derived a scaling expression-a new phenomenologic law, Eq. ( 18), for the flow rate dependency of relative permeabilities.Fractional flow analysis is complementary to energy efficiency analysis.The scope is to provide a clearer image of the combined effects of viscosity, interfacial tension, wettability, pore network geometry, flow intensity, etc., on the structure of interstitial flows and relative permeabilities.The practical implementations relay to the use of a true-to-mechanism phenomenological law for two-phase flow in porous media.Potential applications could be integration in field-scale simulators or machine learning interventions to improve specificity/accuracy/ training; also to improve existing SCAL protocols and associated workload in flow system-or rock-typing applications.
The paper is organized in 6 sections, including this introductory section.Section 2 provides a basic analysis of steady-state, two-phase flow in porous media, and a discussion on the concept of fully developed flow conditions.Section 3 presents the results of extended simulations with the DeProF model algorithm.Reduced pressure gradient values, x, are presented in plane (2D) diagrams; inherent trends are revealed and identified.Section 4, the core section of this work, presents how indicative, generic functional forms of Ca and r, describing the flow rate dependency of the relative permeabilities, are extracted from available data.The concept of intrinsic dynamic capillary pressure is also introduced here.Section 5 derives the corresponding functional form of the energy utilization index, f EU , and the associated locus of critical flow conditions; characteristic invariant properties of the process are also revealed.Discussion, conclusions and suggestions for future research are presented in Sect.6.

Basic Considerations
Consider the simultaneous, one-dimensional concurrent flow of a non-wetting phase and a wetting phase, at flow rates equal to qn and qw , respectively, through a stream tube of fixed cross-section within a porous medium, Fig. 1.The cross-section of the cylindrical control volume may have any arbitrary shape, so long as, it is large enough to be considered homogeneous and side boundaries have no significant effect on the flow structure.
As a result of the externally imposed flow rates, corresponding pressure gradients, (Δp∕Δz) i , i = n, w, are induced upon the two phases.This type of concurrent immiscible flow is conventionally described by the set of phenomenological fractional flow Darcy relations, where indices i = n, w indicate the wetting phase (w) and the non-wetting phase (n), both incompressible, qi are the corresponding flowrates across a porous medium surface of area Ã , perpendicular to the direction of macroscopic flow, Ũi are the corresponding superficial velocities, or flow rate intensities, (Δp∕Δz) i the corresponding pressure differences over any length along the macroscopic flow direction, or pressure gradients along the stream tube, k is the absolute permeability of the pore network and, i are the dynamic viscosities of the two fluids.
Because of bulk viscosity and capillary hysteresis resistances, energy dissipation is manifested as pressure drop along the flow direction.A relative measure of the viscous forces over the capillary forces is provided by the capillary number.Here, we will use the conventional expression for evaluating the capillary number, Ca, as where γnw is the interfacial tension between the two phases.
Please note that the expression, Eq. ( 2), is conventional.It does not take into account the contribution of the viscous forces induced by the flow of both phases, neither those induced by capillarity and wettability, associated with the contact angles of the menisci.The latter, together with the asymmetries and micro heterogeneity in pore geometry, and other intrinsic, hysteresis effects associated with fluctuation dissipation, are the main sources of what is referenced as "capillary resistances".These contributions have a significant impact on the resulting macroscopic pressure gradient(s).In that context, Oughanem et al. (2015) have extended the definition of the capillary number to take into account the effect of wettability.Valavanides (2018b) has proposed to include implicitly the relative magnitude of viscous to capillary forces by considering an effective, appropriately scaled, definition of the capillary number.Yet, that requires an assessment of the flow system over a broad domain of flow conditions.
The flow rate ratio, r, is defined as and is essentially a kinematic property of the flow.
Ca and r are considered as the operational or flow parameters of the process.The set of flow parameters is complemented by the set of system parameters, namely the absolute permeability of the pore network, k , the bulk viscosity ratio, , defined as the physicochemical characteristics of the two phases, e.g., the interfacial tension, γnw , the dynamic, receding and advancing contact angles of the menisci, R and A , and relevant wettability properties, as well as other geometrical, topological and structural characteristics of the pore network.
So long as the process is regulated by imposing a set of flow rate values, qn , qw , across a porous medium surface, Ã , the superficial velocities, Ũn , Ũw , can be con- sidered to be the independent variables of the process.Therefore, without any loss of generality, Ca and r, from Eqs. ( 2) and (3), can be considered as the independent variables in the set of the phenomenological (cause-effect) fractional flow Darcy relations that describe the steady-state process.In that context, the pressure gradients, (Δp∕Δz) i , i = n, w , are the dependent variables.The absolute permeability of the porous medium and the fluid viscosities are considered to be fixed; therefore, the relative permeabilities need to be determined for different flow conditions to form a definite set of cause-effect Darcy fractional flow relations.In the reservoir engineering (2) Ca = μw Ũw γnw (3) r = qn ∕q w = Ũn ∕ Ũw (4)  = μn ∕ μw domain, the applicable standard is to use one of the fractional flows, be it the nonwetting phase, f n , or the wetting phase, f w , correspondingly defined as In the present work, we use as independent variable, the flow rate ratio, r, instead of the fractional flow, because it has the advantage of a more convenient description of the sought physical process, especially in detecting the critical flow conditions of maximum energy efficiency (see Sect. 5 and Valavanides et al. 2016).Switching between flowrate ratio and fractional flow is directly applicable from Eqs. ( 5).
The mobility ratio, , is defined as the ratio of the mobilities, and is essentially a kinetic property of the flow as it relays to the pressure gradients, (Δp∕Δz) i , i = n, w.
When steady-state, fully developed flow conditions are established (see Appendix I), the pressure gradient is common in both fluids and the phenomenological (causeeffect) fractional flow Darcy relations become, Consequently, from Eq. ( 6), it turns out that a basic flow characteristic of steadystate fully developed flow is that the flow rate ratio becomes equal to the mobility ratio, The equivalence between flow rate ratio and mobility ratio in steady-state fully developed flow conditions is a critical characteristic property for the analysis of twophase flows in pore networks.The assumption in Eq. ( 8) was referenced in earlier works (Avraam and Payatakes 1999;Valavanides and Payatakes 2001).Equality of pressure gradients is a trustworthy assumption, provided the flow is fully developed, and saturation gradients are negligible along any REV stream tube.This is equivalent to negligible end-effects if core-scale measurements are referred.On the same line, Standnes et al. (2017) have proposed a mechanistic model for estimating generalized relative permeabilities, whereby they have considered that equal pressure gradients settle in both phases within the representative volume element, neglecting end-effects.The kinematic-kinetic equivalence described by the expression in Eq. ( 8) can also be used in practice as a criterion on whether a specimen of pore network or porous medium can be considered as REV or not, and if fully developed flow conditions are established and to what degree, especially when end-effects are expected to be critical (Karadimitriou et al. 2023). (5)

The Reduced Pressure Gradient for Fully Developed Two-Phase Flows in Porous Media
The conventional fractional Darcy description of immiscible, fully developed, two-phase flow in porous media is provided by Eqs. ( 7).We may recast these equations in a dimensionless form as follows.
We introduce the normative pressure gradient of an equivalent one-phase (saturated) flow of the wetting phase, (Δp∕Δz) 1Φ , at superficial velocity equal to Ũw , We may then use this pressure gradient to normalize the pressure gradient, induced during fully developed two-phase flow in a porous medium, into the reduced (dimensionless) macroscopic pressure gradient, x , given by From the equivalence between flow rate ratio and mobility ratio, Eq. ( 8), at fully developed flow conditions, relative permeability values can be computed directly from the reduced pressure gradient, x, as The above relations comprise the essential definitions of relative permeabilities.When fully saturated flow conditions are met, x = k rw = 1.
Using the reduced pressure gradient, two-phase flow in porous media may be described by a generalized functional form, where A , R , represent the dynamic advancing (A) and receding (R) contact angles of the o/w interface on the pore walls, and x pm is a (sub-) set of (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, chamberto-throat size correlation factors).
The set of physical variables in the parenthesis comprise the set of independent variables.These are classified into two groups.The group of the flow parameters (Ca, r) can be regulated during the process, and the group of the system parameters , A , R , x pm are considered to have fixed values per flow system.The two groups are separated in Eq. ( 12) by a semicolon.The benefit of describing the sought process at fully developed flow conditions in the general form of Eq. ( 12) is the reduction in redundancies associated with the many physicochemical parameters of the flow system, without compromising the consistency of the underlying, true-to-mechanism, flow model.
The functional form in Eq. ( 12) is unknown.The objective of the present work is to extract a generalized, universal expression that can describe the physics of the flow consistently and that can be appropriately adjusted to match the behavior of different flow systems.In the next sections, we will see how an explicit function of the generalized, reduced pressure gradient, Eq. ( 12), has been revealed by matching extensive, systematic simulations, implementing a hybrid mechanistic-stochastic model (the DeProF model), over a broad domain (Ca, r) of fully developed flow conditions for a set of physicochemical flow system parameters , A , R , x pm .Implementing a similar methodology, this functional form can be adapted to describe flows in different two-fluid/porous medium systems.

Prediction of Relative Permeabilities for Steady-State, Fully
Developed, Two-Phase Flow in Model Pore Networks

DeProF Model for Two-Phase Flow in Pore Networks
The DeProF model is a hybrid mechanistic-stochastic model developed to describe immiscible, steady-state two-phase flow in pore networks (Valavanides 2018a).The model is based on the concept of decomposition of the macroscopic flow into prototype flows, hence the acronym DeProF.It takes into account the pore scale mechanisms and the sources of nonlinearity caused by the motion of interfaces, as well as other complex, network-wide cooperative effects, to estimate the conductivity of different classes of model pore unit cells.The different classes of unit cells are partitioned not only on geometric configurations but also on flow configurations, e.g., conducting/non-conducting, fully or partly saturated and containing one or several n/w menisci.Some critical characteristics of the DeProF model comprise: (a) Decomposition of the total flow into three constituent flows, namely Drop Traffic Flow, Ganglion Dynamics and Connected Pathway Flow, each with distinctive sizing of the disconnected fluidic elements of the non-wetting phase.All three comprise an ensemble of interacting fractional flows wandering within a phase space of constituent flows, see points (e) and (f).(b) Dynamic contact angles for the receding and advancing interfaces associated with drainage and imbibition events within pore unit cells.(c) Intermittent flow of fluidic elements is taken into account as a process stochastic property, based on calculated stranding and mobilization probabilities.(d) Scaling-up implemented by applying effective medium theory with appropriate expressions for pore-to-macro-scale consistency for mass and momentum balances over the wetting and non-wetting phases, result in an implicit expression in the form of Eq. ( 12).(e) Determination of the ensemble of physically admissible interstitial flow configurations or system flow microstates, consistent with the externally imposed flow conditions, satisfying Eq. ( 12).(f) Implementation of ergodicity to calculate the expected macroscopic, average reduced pressure gradient, from the ensemble of physically admissible flow microstates, see (e), for different steady-state, fully developed flow conditions and flow systems.
Using the DeProF model, we may obtain a semi-analytical solution to the problem of steady-state, fully developed, two-phase flow in pore networks in the form of Eq. ( 12).The model also predicts key physical quantities that describe the interstitial, bulk and interface structure of the flow: ganglion size and ganglion velocity distribution, fractions of mobilized/stranded populations, specific (per unit volume of porous medium) surface area, velocity and volume fractions of mobilized and stranded interfaces, degree of fragmentation of the non-wetting phase.If the geometry of the pore unit cells is explicitly provided in analytic expressions, the DeProF algorithm is relatively fast, with typical times between 10 and 60 min per (Ca, r)-run on a desktop computer, depending on the userselected flow system properties and flow conditions.

DeProF Model Simulations
Extensive simulations implementing the DeProF algorithm have been carried out over a model pore network with characteristic properties similar to a typical Berea sandstone.Output data were used to derive maps describing the dependence of the flow structure on the independent flow variables, namely the capillary number, Ca, and the flow rate ratio, r (Valavanides 2018a).The simulations span 5 orders of magnitude in Ca (−9 ≤ log Ca ≤ −4) and in r (−2 ≤ log r ≤ 2) .Fluid systems with various-favorable and unfavorable-vis- cosity ratios,  = μn ∕ μw ∈ {0.2, 0.33, 0.67, 1.0, 1.5, 3.0, 9.0, 20} , have been examined.The model flow systems differ only on the value of the viscosity ratio.Properties pertaining to capillarity, i.e., the interfacial tension, wettability, expressed by the dynamic advancing and receding contact angles, as well as the model pore network geometric/topologic characteristics, were kept fixed.The afore mentioned domain of flow conditions has been scanned in 41 steps over Ca, i.e., for values Ca mj = m × 10 j , m ∈ {1, 2, … , 10}, j ∈ {−8, … , −5} , and 21 steps over logr, i.e., for values log Typical predictions of the reduced pressure gradient, x, for different values of Ca and r, and for different N/W viscosity ratio values, , are presented in Fig. 2. Diagrams are shaped into two classes, "iso-Ca" iso-capillary (left) and "iso-r" iso-flow rate ratio (right).
In the following, we will identify and reveal some invariant characteristic properties of the process.

DeProF Model Simulations: Deciphering the Results
In single-phase flows in pore networks, the macroscale pressure gradient scales linearly with the superficial fluid velocity-as long as the flow is maintained below the Forscheimer regime.Flow linearity is also observed for certain two-phase flow conditions, especially so when the momentum and/or energy exchange between the two phases and through n/w interfaces (menisci) is minimum.Such conditions are met when disconnected fluidic elements of the non-wetting phase are stranded and flow of the non-wetting phase is through continuous pathways (fingers), and/or any capillary interactions are negligible when compared to the viscous resistances of the bulk phase flows.
This type of linear behavior is also predicted in the DeProF model simulations for certain flow conditions.In particular, we may note on the x i , log r i iso-Ca diagrams in Fig. 2 that, at the high-end of the flow rate ratio domain, all curves pertaining to constant-Ca values tend to converge asymptotically to an inclined straight dashed line; its functional form is given by the expression Equation ( 13) states that at sufficiently large values of the flow rate ratio, r → ∞, the reduced pressure gradient, x, becomes a linear function of r, by a linearity constant equal to the viscosity ratio.Essentially, Eq. ( 13) describes the decoupling between the wetting and nonwetting phase flows.It can be recovered analytically by writing the definition of the reduced pressure gradient, x, as the ratio of the pressure gradient of the two-phase flow normalized with the pressure gradient that would have been built considering the flow of the wetting phase with superficial velocity equal to the superficial velocity during the two-phase flow, Eq. ( 10).Now, considering the cases of very large values of the flow rate ratio, r → ∞ , the (com- mon) pressure gradient (at steady-state fully developed flow conditions) approaches the following limit (of virtually saturated Darcy flow of the non-wetting phase), Then, at this limit, simply by inserting Eq. ( 14) into the definition of the reduced pressure gradient, Eq. ( 10), we can recover a linear relation between the reduced pressure gradient, x, and the flow rate ratio, r, The higher the values of the capillary number, the earlier (in terms of flowrate ratio values) the approach to that decoupled-flow behavior, and the earlier Eq.( 15) is recovered.Moreover, the logarithmic relation in Eq. ( 13) indicates that the gradient of the asymptote in the (log x, log Ca) diagrams in Fig. 2 should be 1:1, irrespective of the system properties.This can be verified by observation across all the x i , log r i iso-Ca diagrams, irrespective of the value of the viscosity ratio.
Considering the x i , log Ca i iso-r diagrams, aligned on the right columns in Fig. 2, we observe a similarly interesting trend.At the low-end of the capillary number domain, all iso-r curves tend to converge asymptotically to a straight line, illustrated by the dashed thick line with (negative) inclination.The inclination of the asymptote is −1 ∶ 1 , across all the x i , log Ca i iso-r diagrams, irrespective of the value of the viscosity ratio, .Again, decoupling between the flows of the wetting and non-wetting phases occurs at very low-Ca values.The lower the value of flowrate ratio, this asymptotic behavior is observed from higher Ca values.Moreover, we may verify by direct measurements that, in all diagrams, irrespective of the value of , the low-Ca asymptote crosses the x = 1 line (pertaining to "virtually saturated flow") at the same value of the capillary number, log Ca pm = −4.75 .That is indicated by a large red cross ( ).In that context, the functional form of the low-Ca asymptote is given by Also note, in Eq. ( 16), that the value of the constant, Ca pm , across all the fluid systems examined, is Therefore, this value of Ca pm is considered to be a characteristic value (a property) of the fluid system.It depends on capillarity related characteristics, i.e., the tripartite interaction between interfacial tension, wettability (dynamic advancing and receding contact angles) and the pore network geometry.Should we have had the simulations rerun, with either a different wettability or different pore network geometry, that value would had changed.The open question is how sensitive the behavior of the fluid system/process would be to such changes, not only in simulations of virtual systems but also in applications in natural systems.
We see that, at the extremities of the flow conditions, the conventional assumption of saturation dependency of relative permeabilities is recovered.Yet, this is only a small part of the entire domain of flow conditions.For the rest of the flow domain, the broad, central region of moderate Ca values, where viscosity-induced resistances compete with capillarity-induced resistances, relative permeabilities are strongly nonlinear functions of the flow intensities of the wetting and non-wetting phases.In the following section, we will see how a universal functional form, x(Ca, r) , can be derived to fit the data, x ij Ca i , r j , predicted from the DeProF model simulations.
4 The Scaling Function for the Flow Rate-Dependent, Reduced Pressure Gradient

Introduction of the Scaling Function Model
We will continue with deriving a functional form for the reduced pressure gradient, x(Ca, r) , in order to best fit the data, x ij Ca i , r j , produced from the DeProF model simulations for 7 different systems of n/w viscosity ratio values, ∈ {0.33, 0.67, 1.0, 1.5, 3.0, 9.0, 20} .The diagrams in Fig. 2 are representative data sets x ij , for various viscosity ratio values.
We may start by examining the self-similar structure of the x ij values on the iso-Ca, diagrams stacked on the left column diagrams, in Fig. 2. In parallel we will observe the inherent asymptotic trends described by Eqs. ( 15) and ( 16).
We may observe a systematic trend, repeating for every fixed value of Ca.As log r is progressively increased, x starts from a constant value and then it progressively increases-following a power law trend-toward attaining a linear dependence on r.This trend can be described with the following functional form In Eq. ( 18), the constant value at the low-r regime, when r → 0 , is represented by the term A, while the linear trend toward the asymptote line at the high-r regime, by r .We may also observe on the left column diagrams in Fig. 2 that, the constant value of x attained at the low-r regime decreases with increasing Ca.Therefore A is correctly considered to be a function of the capillary number, A = A(log Ca) .Moreover, all the iso-Ca lines converge asymptotically as r → 0 , at the same line, x = r .At those condi- tions, any capillarity effects, expressed by A(log Ca) , should decrease with increasing Ca.This indicates that A(log Ca) should be a strongly decreasing function of Ca.The same self-similarity repeats consistently for different viscosity ratio values, .Focusing now on the right-column diagrams in Fig. 2, we may observe the trend of the iso-r lines into converging asymptotically, as Ca → ∞ , to the horizontal axis at value x = 1 , indicating reduced pressure gradient values pertaining to saturated flow conditions.That is consistent with the definition of the reduced pressure gradient, Eq. ( 10), and gives an asymptotic behavior described by, When flow conditions shift from intermediate to relatively high-Ca values, the interstitial flow progressively mutates from prevailing ganglion dynamics to a flow mixture rich in drop traffic flow and connected pathway flow.With further increase in Ca the non-wetting phase flow becomes bimodal, comprising connected pathways and drop traffic flow, whereby the relatively large drops of the non-wetting phase shutter into smaller droplets.The latter progressively loose contact with the pore walls, and the disconnected flow mutates to an emulsion-type flow mode.Consequently, in Eq. ( 19), the effects of menisci and capillary resistances fade away, A → 1 , and viscosity becomes the regulating factor, r .Details on the interstitial flow structure mutations as flow condi- tions intensify can be found in Valavanides (2018a).
The functional form, A(log Ca) , describing the effect of the n/w menisci at different Ca values is revealed in the following paragraph by simple, least squares fitting.

Recovery of the Kernel Function A log Ca
We will see how the kernel function A(log Ca) in Eq. ( 18) can be recovered for the particu- lar flow system we examined.In principle, we need to have a sufficient amount of relative permeability values extending across a broad domain of flow conditions Ca and r.In the present paradigm, we will tap on the extensive DeProF model simulations (virtual experimens), covering a broad domain of flow conditions, log Ca and log r , and for a variety of viscosity ratio values.Wettability and pore network geometry remained fixed across all simulations.
The reader who is not interested in following the technicalities of the fitting methodology may skip the details and directly refer to the result expressed in Eq. ( 26).
The functional form of the kernel function, A(log Ca) , depends on the particular fluid system examined.Similar paradigms, in extracting A(log Ca) from laboratory experiments, are provided in Valavanides et al. (2020) and Karadimitriou et al. (2023).

Extraction of the Kernel Function, A log Ca , by Fitting
We will implement least squares fitting to recover the A(log Ca) kernel function in Eq. ( 18).For any data set pertaining to a fixed-Ca simulation, log Ca i , we detect the value A i = A log Ca i , that minimizes the sum of the squares of the fitting errors in each of the iso-Ca data sets (each set depicted with connected line in the diagrams exhibited on the left columns in Fig. 2).This is expressed as where N i is the number of detected points (x ij ) pertaining to the constant capillary number value, Ca = Ca i , for N i different flowrate ratio values, r j , ( j = 1, 2, … , N i ≤ 21 ).Please note that the last expression for A i in Eq. ( 20) is the average We recall here that each set A i = A log Ca i , 1 ≤ i ≤ 41 pertains to a particular value of the n/w viscosity ratio, .So, for each examined value of the viscosity ratio, we get a different set comprising 41 pairs of values, A i , log Ca i , 1 ≤ i ≤ 41.
We repeat the procedure for each of the 7 viscosity ratio values, ∈ {0.33, 0.67, 1.0, 1.5, 3.0, 9.0, 20} that we have examined in the DeProF model simu- lations.These values are plotted in the diagram Fig. 3a.The 7 sets of coefficient values, (20) , are depicted using 7 types of markers.Each type of marker indicates 41 data values A i , 1 ≤ i ≤ 41 , pertaining to the data fitting, Eq. ( 20), for the correspond- ing constant-Ca i runs.It is obvious that the 7 sets of A i , log Ca i , 1 ≤ i ≤ 41 points are lined-up on a tight formation.We may assess how the "tightness" is distributed over the 41 log Ca i values by calculating the corresponding variances.
The 41 Ca i variances over the 7 viscosity ratios are evaluated by The 41 variances per Ca i of each set of log A values, pertaining to the 7 different vis- cosity ratios, are presented in the diagram, Fig. 3b.Variances pertaining to the capillarity dominated regimes, log Ca < −6, are a bit higher than those in the capillary/viscous flow regimes, −6 < log Ca .Nevertheless, as the variances are quite small, we may calculate the average over the viscosity ratio values, ⟨A i ⟩ m , as The 41 averages calculated from Eq. ( 23) are plotted on the diagram Fig. 3c with 41 red dot markers ( ).
On the diagrams Fig. 3a and b, we may observe two distinct regions along the log Ca domain.At the low-Ca regime, the markers are almost aligned on a straight line with negative inclination, to be calculated by fitting.As the capillary number is increased, the ( 22) .33, 0.67, 1.0, 1.5, 3.0, 9.0, 20} 〈log 〉 , Eq. ( 23) log (log ), Eq. ( 26 2 , Eq. ( 22) 2 , Eq.( 22) These are fitted by a function A(log Ca) (black, dash-dot curve), of the form Eq. ( 26) markers tend to progressively attain the lowest physically admissible value (log A → 0) , and in doing so, they align asymptotically to the horizontal direction.Focusing on the diagram Fig. 3c, each of the 1 ≤ i ≤ 41 red dot markers ( ) indicates the corresponding average ⟨A i ⟩ m over the 1 ≤ m ≤ 7 viscosity ratios, calculated from Eq. ( 23).The red dot markers ( ) align along the inclined continuous line in black; the latter meets the horizontal axis at an obtuse angle at point C, at abscissa log Ca 0 .
We may then fit these markers by a bilinear function of the form where C 0 is a measure of the approach to the two asymptotes (the higher this value, the the distance from the asymptotes); −1∕C 3 is the gradient of the inclined asymp- tote; log Ca 0 is the abscissa of the intersection of the two asymptotes (at point C).
Best fitting of the ensemble of log A values, for the entire set of viscosity ratios exam- ined in the simulations, with Eq. ( 24), gives the values of the coefficients presented in Table 1.The curve fitting is plotted by the black dash-dot line in Fig. 3c.
The red dashed line, pivoted at log Ca pm = −4.75 , ( ) with negative gradient 1:1, see Eq. ( 17), is the asymptote delineated in the right column diagrams in Fig. 2. That is an indication of the self-similarity properties of the kernel function A(log Ca).
Returning to Eq. ( 24), log A is an implicit function of log Ca .We need to solve the corresponding cubic equation in terms of log A, The solution to the previous cubic equation is derived in "Appendix II" and is given by the following expression, where One may select to choose between the implicit or explicit expressions for log A(log Ca) , Eqs. ( 24) and ( 26), respectively.Note, the expression provided in Eq. ( 26) pertains to the particular fluid model system examined in the simulations and has no universal application.Evidently, each flow system, comprising the two fluids and the ( 24) Table 1 Estimated values of the A(log Ca) function coefficients, Eq. ( 24) porous medium, has its own "finger-print" expression for log A(log Ca) .Similar para- digms in extracting A(log Ca) from laboratory measurements of different flow systems can be found in the literature (Valavanides et al. 2020;Karadimitriou et al. 2023).In the next section, we will interpret the physics underlying the A(log Ca) curve.

Physical Interpretation of the Kernel,
A log Ca : The Intrinsic Dynamic Capillary Pressure Function As described in the previous paragraph, by fitting the data obtained from the simulations, Fig. 3, we get the functional form of A(log Ca) .Obviously the values of the constants, Ca 0 , C 0 and C 3 , do not change with the viscosity ratio, at least not substantially and to the extent that algorithmic deficiencies of the DeProF model are unavoidable.Therefore, simulations indicate that A(log Ca) is invariant to changes in viscosity ratio.Now, rearranging Eq. ( 18), we get In that format, the function A(log Ca) is describing the effect of the n/w menisci on the transport configuration of the flow system at different flow conditions, i.e., at different flow intensities, or different Ca and r values.In Eq. ( 28), by subtracting the effects of flow resistances due to bulk viscosities, r , from the reduced pressure gradient, x, we get the relative contribution of the dynamic capillary resistances to the flow.Therefore, the values of the constants, Ca 0 , C 0 , and C 3 , in Eq. ( 18) are expected to change with wettability (dynamic advancing and receding contact angles) and/or pore network geometry.
The dimensionless function, A(log Ca) , may be given the name Intrinsic Dynamic Cap- illary Pressure-IDCP function.As such, the IDCP function is a characteristic property of the flow system, comprising the two fluids and the porous medium.It depends on capillarity related characteristics, i.e., the tripartite interaction between interfacial tension, wettability (dynamic advancing and receding contact angles) and the pore network geometry.The term dynamic in IDCP is used to indicate that capillary pressure is evaluated at different flow conditions, not in static equilibrium.
The IDCP curve plotted in the diagram in Fig. 3c (black dash-dot curve) represents the particular fluid system examined in the simulations (Eqs.( 26), ( 27) and Table 1).Should we have had the simulations rerun-with either different wettability properties or different pore network geometry-its functional form, Eq. ( 28), would have changed.Any changes would be due to changes in the constituent physical variables of Eq. ( 28) and, most critically, to changes in the structure of the response of the system, x(Ca, r) .That would be the case if, the absolute permeability of the network in saturated flow conditions remains fixed, and changes would only pertain to wettability aspects.Albeit the IDCP curve is revealed from REV-virtually ex-core-measurements, it contains the essential information on the interaction between the structure of the network (associated with capillarity) and the motion of the n/w menisci.In that context, the IDCP curve surfaces-up, to a superficial, excore observer, the latent information associated with the ensemble of flow mechanisms at the pore scale, and the corresponding interstitial flow structures.
In the current paradigm of the IDCP curve in diagram Fig. 3c, we may identify 3 distinct regions.Starting from the capillary flow domain at relatively low-Ca values (log Ca < −5) , the IDCP curve is decreasing in an almost linear fashion, indicating a progressive mobilization of stranded fluidic elements, mainly the larger size ganglia.Then, there is a transition region of mixed capillary-viscous flow at moderate-Ca values (−5 < log Ca < −4) , where (28) A(log Ca) = x(Ca, r) − r the IDCP curve gradient progressively decreases.This is because the remaining stranded ganglia of smaller size has a smaller probability of mobilization.Then, moving to the viscous dominated flow regime of high-Ca values (−4 < log Ca) , the transitory behavior of the flow system turns into a smooth asymptotic behavior toward saturated flow (A → 1) , as explained in the paragraph following Eq.( 19).18) and ( 26) with parameter values from Table 1.The fitting is quite consistent and accurate, both qualitatively and quantitatively

Maps of Reduced Pressure Gradient and Relative Permeabilities
Using the scaling function Eqs. ( 18) and ( 26) with coefficient values from Table 1, we may reconstruct 2D, iso-Ca and iso-r diagrams, similar to those presented in Fig. 2. The diagrams in Fig. 4  Diagrams of x(Ca, r) and combined diagrams of k rn (Ca, r) and k rw (Ca, r) for different flu- ids with viscosity ratio values, ∈ {0.33, 0.67, 1.0, 1.5, 3.0, 9.0, 20} , are presented in Figs. 5  and 6, respectively.The diagrams in Fig. 4, columns (b) and (d), are formed by intersections of the x(Ca, r) 3D surface with iso-Ca and iso-r planes, perpendicular to the log Ca and the log r axes, respectively.These are merged to draw the 3D diagrams in Fig. 5.The combined relative permeability diagrams of k rn and k rw in Fig. 6 are drawn using Eq. ( 11).
We may focus on the intersection of the relative permeability surfaces, forming an S-type plane curve, the locus of cross-over relative permeability values, k rx = k rn = k rw .We may observe its projections on the 3 mutually orthogonal planes of reference, k rn = k rw = 0 , log Ca = −8 and log r = −2 .The two projections on the planes k rn = k rw = 0 and log Ca = −8 are aligned on the two straight lines, r = r x = −1 , perpendicular to each other.That is a consequence of the equivalence between the flow rate ratio and the mobility ratio, Eq. ( 8), in fully developed flow: At the cross-over relative permeability flow conditions, the flow rate ratio and the mobility ratio become reciprocal to the viscosity ratio, The projection of the cross-over relative permeability curve on any fixed-Ca plane is invariant to viscosity ratio changes, so long as both wettability and pore network geometry remain fixed.Therefore, cross-over relative permeability conditions, k rn = k rw = k rx at r x = 1∕ , will remain unaffected when flow systems with different viscosity ratio but with the same wettability and pore network geometry are examined.
This latent (or less-expected) characteristic property, comes as a direct consequence of the invariance of the IDCP curve, A(Ca) , to viscosity ratio changes (see Fig. 3a), and altogether may lead to the conjecture that the functional dependence of the cross-over relative permeability on Ca is a pure capillarity dominated system characteristic.As such, it could be used as an effective, pore network structure identity.Because of its potential implications in improving SCAL technology, this conjecture, merits systematic laboratory verification.

Energy Efficiency and Critical Flow Conditions
The energy efficiency or energy utilization index, f EU , provides a reduced measure of the energy efficiency of the process.It is defined as the flow rate of the non-wetting phase per unit of the total hydraulic power spent, or equivalently, provided externally to the flow system by the driving/injection pumps to maintain two-phase flow at set flow conditions, Ca and r (Valavanides et al. 2016).It is used as a flow analysis and system characterization variable.The energy efficiency index may be evaluated in terms of macroscopic (REV) measurements, (29) ∀r = r x ∶ k rn Ca, r x = k rw Ca, r x Eq. (8) ⇒ r x = 1∕ Fig. 5 3D plots of the reduced pressure gradient, x, as a function of the capillary number, Ca, and of the flow rate ratio, r, based on Eqs. ( 18) and ( 26) with parameter values from Table 1, and for and for various n/w viscosity ratio values, ∈ {0.33, 1.0, 1.5, 20} Fig. 6 3D plots of the relative permeability for the wetting and non-wetting phase, k rw (blue) and k rn (light red), as a function of the flow rate ratio, r, and the capillary number, Ca, based on Eqs. ( 11), ( 18) and ( 26), with parameter values from Table 1, and for pairs of fluids with different n/w viscosity ratio values, ∈ {0.33, 1.0, 1.5, 20} .The plane, S-shaped curve, formed by the intersection of the two relative permeability surfaces, the locus of cross-over relative permeabilities is invariant to viscosity ratio changes The energy efficiency index has an important characteristic.For every fixed value of the capillary number, Ca, there exists a single value of the flow rate ratio, r*, for which the energy efficiency index, f * EU , attains a maximum value.This maximum is local; changing Ca to Ca′, the new maximum is met at another unique value r * � .Therefore, for every flow system, comprising the two fluids and the porous medium, a unique locus of energy efficiency maxima is formed, r * (Ca) .Flow conditions matching the r * (Ca) locus are called critical flow conditions (CFC) (Valavanides 2018b).Critical flow conditions can be revealed by special core analysis (Valavanides et al. 2016;2020;Karadimitriou et al. 2023) provided that measurements are taken over an extended range of fully developed flow conditions.
Given the scaling function for the reduced pressure gradient, x(Ca, r) , Eq. ( 18), we may obtain, through Eq. ( 30), the analytic expression for the energy efficiency index, f EU (Ca, r) .We may also derive the analytic expression for the locus of CFC by its con- ceptual definition, whereby importing the expression for the x(Ca, r) scaling, Eq. ( 18), we get the CFC locus.
The locus of critical flow conditions is configured by two separate factors, the viscosity ratio, , incorporating the bulk viscosity effects, and the IDCP function, A(Ca) , incorporating the capillarity effects.That is something expected as f * EU (Ca, r) is essen- tially the minimum energy partitioning between viscosity-and capillarity-induced dissipation.
Plots of the energy efficiency index and critical flow conditions as a function of Ca and r, Eqs. ( 30), ( 18) and ( 26), for various viscosity ratio systems, are presented in Fig. 7.In every diagram, the 3D surface with the iso-f EU contour lines indicates the dependence of the energy efficiency index on log Ca and log r .The thick red line indicates the trajec- tory of maximum energy efficiency values, f * EU (Ca, r) = f EU (r * (Ca)) , based on Eq. ( 32).On that trajectory, magenta colored ball markers are used to pinpoint the trajectory of f * EU (Ca, r) over order of magnitude shifts in Ca, [ log Ca ∈ {−4, −5, −6, −7, −8}].
The projection of the trajectory of maximum energy efficiency, f * EU (Ca, r) , showing the locus of critical flow conditions, r * (Ca) , on the log Ca × log r plain, lays-barely visible-under the f EU (Ca, r) surface.To show the locus clearly, the f * EU (Ca, r) curve is re-projected, on a horizontal plane (Ca × r) located just above the 3D surface.Very thin red lines connecting the f EU (Ca, r) surface ridge to its new projection appear as a cylin- drical curtain.The new projection depicted by a red dash-dotted curve lays at a height equal to the nominal value of the maximum attainable energy efficiency-the "ceiling of efficiency"-for the particular fluid system-better, for the particular value of the viscosity ratio (Valavanides 2018b).That value corresponds to the critical flow rate ratio, specifically, r * ∞ = 1∕ √ , and is theoretically estimated as The straight line asymptote, r = r * ∞ = 1∕ √ , is depicted with dashed red line and lays, barely visible, under the f EU (Ca, r) surface.Both of the aforementioned nominal values, r * and f * EU∞ , are estimated by considering the asymptotic limit at the high-end of the capillary number domain, whereby the flow resistance due to the effects of the n/w menisci and phenomena associated with capillarity, become negligible when compared to bulk viscosity flow resistances (Valavanides 2018b, Appendix).
The shapes of the energy efficiency surface ridge, f EU (r * (Ca)) and of the locus of criti- cal flow conditions, r * (Ca) , become more clear if we observe their projections on the three reference planes, Fig. 7.In particular, the projection of the maximum energy efficiency on the log r = −2 plane attains an S-shaped curve, confined asymptotically between the two horizontal lines, f EU = 0 and . Similarly, the projection on the log Ca = −8 plane attains a smooth, half-U shaped, bilinear curve, delimited by two mutu- ally perpendicular asymptotes, f EU = 0 and log r = log r * ∞ = −0.5 log .Α set of three or four parameters may be used to describe, at an adequate level of approximation the particular S-shape projection.Furthermore, we may normalize the S-shape expression-associated with any particular two-fluid system-with the corresponding value of the maximum attainable efficiency (a n/w fluid system property) and, in doing so, we may also scale the Ca values, conventionally defined by Eq. ( 2), to a universal, system-reduced capillary number, defined as Ca R = Ca C 1∕z , with parameters C and z infusing the combined effects of wettability and viscosity ratio of the two fluid/porous medium system (see Valavanides 2018b, Eq. ( 54)).
In such case, all the S-shape curves, pertaining to the same porous medium but with different viscosity ratios, are expected to collapse into a single, normalized S-shape, identical to the S-shape pertaining to the system with equal viscosities, = 1 , as described in the same paper.This way the combined effect of the two bulk viscosities (normalized by the effect of the viscosity ratio) is decoupled from the combined effects of wettability and pore network geometry.The resulting normalized expression of the S-shape function could serve as the "phenotype" of the pore network for a particular wettability.
Τhe previous diagrams may be integrated into a universal relative permeability and energy efficiency diagram for fully developed flow in pore networks.Two diagrams, pertaining to indicative cases, with viscosity ratio values, = 0.67 and = 1.5 , are presented in Fig. 8.The diagrams cover a domain of flow conditions for −8 ≤ log Ca ≤ −3.6 and −2 ≤ log r ≤ 2 .Note, the diagrams are drawn from analytical expressions derived in the previous sections.
Plots of the relative permeability values as a function of log r are presented as tra- jectory curves for 5 fixed-Ca values, log a i ∈ {−8, −7, −6, −5, −4} .These trajectories are the intersections of the 5 fixed-Ca planes with the relative permeability surfaces, k rn (Ca, r) and k rw (Ca, r) , plotted in Fig. 6.The relative permeability plots for the non-wetting phase are drawn in red color and those drawn for the wetting phase in blue color, in consistency with the coloring of the parent surfaces.Markers of different shape/size have been assigned to each fixed-Ca pair of relative permeability curves.The five cross-over relative permeability values at the crossings of the corresponding five pairs of relative permeability curves, are indicated with larger balls (in gray color), while the entire set of cross-over relative permeability values, is depicted by the black dashed, S-shape plane curve.All cross-over points correspond to the same flow rate ratio value, r x ; this is why the S-shape curve lays on the r = r x = 1∕ plane.For the indicative n/w viscosity ratio values, the fixed-r planes correspond to log r x = − log 0.67 = 0.1739 on the top diagram ( = 0.67) and to log r x = − log 1.5 = −0.1760 on the bottom diagram ( = 1.5) .The projections of the cross-over relative permeability surfaces on the two reference planes are depicted with gray, dashed lines; in particular, with a straight, thin gray line on the log Ca = −8 plane and an S-shape, dash gray line on the log r = 2 plane.Plot of the energy efficiency as a function of the flow conditions, f EU (Ca, r) , Eq. ( 30), is shown by the gray color scale 3D surface.Contour lines and mesh lines correspond to the scaling and grid lines on the reference planes.The ridge of the energy efficiency surface, indicating the local maxima attained at critical flow conditions, is delineated by the continuous, thick red curve.The intersections with the fixed-Ca planes, log Ca i ∈ {−8, −7, −6, −5, −4} , are marked with violet ball markers.The projection on the z = 0 reference plane delineates with the red curve the critical flow conditions, r * (Ca) ; it is partly visible as it is obscured by the f EU surface.
Projections on the log Ca = −8 and log r = 2 reference planes are also plotted to depict the inherent configuration characteristics of the maximum energy efficiency curve.The projection on the z = f * EU∞ = 1∕ � 1 + √ � 2 plane, depicted by the red, dash-dotted plane curve provides a clearer layout of the locus of critical flow conditions.Again, as in the diagrams in Fig. 7, the drop lines form a cylindrical curtain, starting from the ridge of the f EU = (Ca, r) surface up to the 'ceiling of efficiency' (at the maximum level of attainable efficiency).The corresponding nominal values are: The energy efficiency map of any two-fluids /porous medium flow system in Fig. 7 or Fig. 8, provides a normative framework for assessing /evaluating the viscous/capillary character of any flow within that system, based on the imposed values of capillary number and flow rate ratio.
This type of diagram is an improved version of the classical phase diagram, proposed by Lenormand et al. (1988) and Lenormand (1990), coarsely depicting the flow structure map in terms of the capillary number and the viscosity ratio, as well as its 3D extension proposed by Ewing and Berkowitz (2001) to account for gravity effects.

Conclusions and Discussion
A generic form of flow rate dependency of the relative permeabilities for steady-state, fully developed two-phase flow in porous media is proposed.The functional form is the result of a two-stage derivation, i.e., systematic simulations and data fitting, extending over a broad range of flow conditions.The latter is expressed in terms of the capillary number, Ca, and the flow rate ratio, r.A variety of two-fluid viscosity ratio systems, , are examined.
In the first stage, virtual experiments of steady-state, fully developed, concurrent flow of two immiscible fluids within a homogeneous model pore network were performed implementing the DeProF model algorithm.The model pore network corresponded to a virtual core of infinite length in the absence of any end-effects.The systematic simulations examined an extended range of flow conditions, spanning 5 orders of magnitude over the capillary number, Ca, and the flow rate ratio, r, and for a variety of viscosity ratio systems, .The main product of the simulations was a dense matrix of reduced pressure gradient values for different values Ca i and r j , x ij Ca i , r j .The reduced pressure gradient measures the relative increase in the pressure gradient of the mixed flow against that of a saturated flow of the non-wetting phase at the same flow rate (superficial velocity).By definition, in fully developed flow conditions, the reduced pressure gradient is equal to the inverse of the relative permeability of the non-wetting phase.The system properties associated with interface/capillary phenomena, i.e., interfacial tension, wettability (dynamic advancing and receding contact angles) and pore network structure, were kept fixed.In that context, it was possible to observe the effects of viscosity, independently from the effects of capillarity.Apart from reduced pressure gradient charts, the systematic simulations produced a complete description of the interstitial flow structure.
In the second stage, combining the essential physics underlying the process with data fitting on the charts produced by the simulations, we derived a universal functional form to describe the dependence of the reduced pressure gradient, x, in terms of the capillary number, Ca, and the flowrate ratio, r.The proposed functional form extends over orders of magnitude in Ca and r, and it is expressed as whereby the effects of capillarity and viscosity are partitioned in two decoupled terms.The linear term describes the fractional contribution of the bulk viscosities through the N/W viscosity ratio, .The nonlinear term, A(log Ca) , accounts the contribution of capillarity at different flow conditions.That nonlinear term is described as the Intrinsic Dynamic Capillary Pressure function.The relative contributions of those two terms depend on the flow conditions.The flow conditions that minimize the total energy dissipation per unit of nonwetting phase flow rate are called critical flow conditions.
The functional form of the flow dependency of relative permeabilities is recovered analytically from the reduced pressure gradient function, by implementing the conventional formulation of the phenomenological fractional Darcy law for fully developed flow, Eqs.(11).The functional form of the energy efficiency and the locus of critical flow conditions can also be derived analytically.
The Intrinsic Dynamic Capillary Pressure curve, formed by the ensemble of reduced, dynamic capillary pressure values at different flow conditions, accounted by A(log Ca) in Eq. ( 34), describes the transition of the capillarity effects on the average flow, between the extremities of the domain of flow conditions, i.e. from capillarity-dominated flows to viscous flows.The IDCP curve is unique for every flow system, and, in that context, it can be used as the identity of the system.Moreover, the IDCP curve is directly related to the physicochemical properties of the two-fluid system and its affinity to the pore network characteristics associated with dynamic capillary phenomena (geometry, topology, wettability).
One of the potential applications of the proposed scaling is in SCAL.That needs to be considered on the basis of established fully developed flow conditions and endeffects.The DeProF model is built around the equality of pressure gradients in both phases.That condition is attained when fully developed flow is in place.That is not a condition always attained in SCAL measurements, as evolution of the interstitial flow depends on the structure of the pore network, the wettability of the fluid system and the flow conditions (the intensity of the co-injection).In general, end-effects are (34) x(Ca, r) = A(log Ca) + r pronounced when the role and effects of capillarity-associated phenomena are not suppressed by viscous resistances.The proposed methodology can be applied in SCAL but the results will be conditionally biased by end-effects.Yet, implementing a few checks, some of these simple, e.g., on the equality between imposed flow rate ratio and measured mobility ratio, and some a bit more complicated, e.g., energy efficiency vs flow conditions, we may detect the presence of end-effects.We have seen that capability in on-going research (Valavanides et al. 2020;Karadimitriou et al. 2023), whereby the effective IDCP functional form has been recovered implementing a similar methodology in fitting the ex-core measurements.Yet, the missing part(s) of the puzzle is on evaluating the extent/severity of end-effects and/or being able to cope with end-effects in unbiasing the ex-core measurements to get a concise and accurate (trustworthy) relative permeability map.
Another interesting application of the proposed approach is in revealing the structure of the interstitial flow.A convincing paradigm on how the IDCP curve may be used as a forensic tool in revealing the structure of the interstitial flow is the systematic SCAL study performed on a core of natural Clasach sandstone (Valavanides et al. 2020).An analysis of the measured relative permeabilities, similar to that described in Sect.4.3 revealed a sharply decreasing intrinsic dynamic capillary pressure function, A(log Ca) , over moder- ate-to/high-Ca flow conditions.That seemed to be contradicting the smooth, asymptotic transition to zero (log A → 0) depicted in Fig. 3a.The sharp decrease in A at moderate/ high-Ca flow conditions can be associated with a sharp reduction in capillary resistances.This could be attributed on progressive eradication of n/w menisci contacting the pore walls.And this could only happen if the small disconnected fluidic elements of the nonwetting phase (ganglia and large drops) break-up into tiny droplets and the immiscible flow of larger drops mutates into some kind of emulsion flow.Reviewing the conditions of the experiments, it turned out that the suspected emulsification was in fact the result of adding surfactant to suppress the interfacial tension in order to impose flow conditions at very high-Ca values.The particular structure of the intrinsic dynamic capillary pressure curve, revealed the onset of emulsification at high-Ca flowrates, an intrinsic flow characteristic, from ex-core measurements.The paradigm indicates the forensic capabilities of the IDCP curve.
Overall, the presented flow analysis is based on a set of appropriately reduced variables, eliminating redundancies, yet without compromising the specificity associated with the description of the process and/or the flow system.The derived analytical expressions of flow rate dependency fit surprising well the predicted values from systematic simulations with the DeProF model.They also provide an additional proof of consistency with respect to the existence of the universal relative permeability and energy efficiency map proposed by Valavanides (2018b), based on an energy analysis of the flow process, from first principles, phenomenological observations and independently of modeling assumptions.
In addition, generalized specific characteristic invariants, pertaining to the two-fluid and porous medium flow system, were revealed: the locus of Critical Flow Conditions, the locus of cross-over relative permeability conditions, as well as the Intrinsic Dynamic Capillary Pressure curve.All of these can be used as flow system identities at different degrees of specificity, i.e., as "genotype" of the two-fluid/porous medium flow system or "phenotype" of the porous medium.This could be a strong tool in developing flow systemor rock-typing methodologies and techniques.
Moreover, the ability to handle a functional form of the relative permeabilities, accounting for a mechanistic, flow rate dependency, may improve the specificity of field-scale simulators, or increase the efficiency in machine learning interventions.
Future research should address how capillarity-related system characteristics, in particular dynamic wettability or network geometry, may affect the loci of Critical Flow Conditions and cross-over relative permeabilities or the Intrinsic Dynamic Capillary Pressure curve.In that context, sensitivity analysis needs to be run on virtual systems, model pore networks or natural cores.Another potential application is in assessing the extent of endeffects in core analysis from ex-core measurements.and consequently a common macroscopic pressure gradient, (Δp∕Δz) , is established.
Establishment of fully developed flow is a critical problem in special core analysis (SCAL) as the measured values of the physical quantities (pressure gradients, saturation and flowrates), that are necessary to evaluate the true phenomenological relation at the REV scale, should not be biased by the, so called, "end-effects".The latter may be substantial or not, depending on the extent of the non-fully developed entry-flow zone, that, on its turn, depends on the physicochemical properties of the two-fluid and porous medium flow system and the imposed flow conditions.In a nutshell, end-effects-if substantialdistort the equivalence between the laboratory examined core flow and the flow within a REV in-situ the formation.In that context, end-effects distort-to a degree-the trueness of the REV constitutive relations, and the specificity of any associated simulations.The effect tends to be significantly reduced, or even negligible, when relatively higher flowrates are imposed, as the non-fully developed length is significantly reduced.In that sense, at higher flow rates, saturation gradients along the core are significantly reduced and saturation attains some uniform value along the core.

II Determination of the Roots of the Cubic Eq. (25)
We will derive the roots of the cubic Eq. ( 25).Without any loss of generality we transform log Ca = x and log r = y , to get a simpler form as where C 3 , C 0 > 0, x 0 < 0, The discriminant of the cubic equation, is positive when and there are 3 distinct real roots (irreducible case).Therefore, we must avoid Cardano's formulae since, in this case, we cannot express the roots in terms of real radicals.For = x 0 , we get For x ≠ x 0 , we implement the transformation, to obtain the cubic equation in reduced form, Since < 0 , we can follow Vieta's method (1591) to solve Eq. ( 12) for x ≠ x 0 .
In general, if we set we obtain three real roots,

Fig. 1 A
Fig. 1 A REV element of one-dimensional flow within a cylindrical, infinitely long, porous medium tube (top).The schematic drawings are successive scale-down depictions of the mixture of the wetting phase (WP, blue) and the non-wetting phase (NWP, red) within the interstitial flow (inset, top), the disconnection of the non-wetting phase (bottom) and the associated hysteresis due to the different contact angles, advancing, bottom left, and receding, bottom right.(The inset is a snapshot from the Avraam and Payatakes (1995) flow experiments in model pore networks, courtesy of D.G. Avraam.)

Fig. 2
Fig. 2 Paired diagrams of predicted reduced pressure gradient values, x, for different values of capillary number, Ca, and flow rate ratio, r.Each pair corresponds to a n/w viscosity ratio value, ∈ {0.33, 1.0, 1.5, 20} .Small markers depict raw values predicted by DeProF simulations (Valavanides 2018a).Markers are connected by straight line segments in groups of iso-Ca values (left parts of diagrams) and iso-r values (right parts of diagrams).In every diagram, the thick dashed line depicts an asymptote, the large blue circle and the red cross indicate characteristic values (see text) For log Ca ≪ 0, log x = log Ca pm − log Ca ⇔ x = Ca pm Ca (17) log Ca pm = −4.749⇔ Ca pm = 1.782 × 10 −5

Fig. 3 a
Fig. 3 a Values of the coefficients A i , 1 ≤ i ≤ 41 , used in Eq. (23) for the 41 constant-Ca runs.Each type of marker indicates a set of 41 data values Ca i , A i and pertains to a system with different viscosity ratio values,  = μn ∕ μw .b Blow-up detail of the calculated variance in each of the 41 log log A i values over the 7 viscosity ratios in (a).(c) The corresponding 41 -averaged values, ⟨A⟩ i , indicated by red dot markers ( ).These are fitted by a function A(log Ca) (black, dash-dot curve), of the form Eq. (26)

Fig. 4
Fig. 4 Reduced pressure gradient values, x, for different values of capillary number, Ca, and flow rate ratio, r, and for various n/w viscosity ratio values, ∈ {0.33, 1.0, 1.5, 20} .Diagrams in columns (a) and (c) plot the simulation results, Fig. 2. Diagrams in columns (b) and (d) show plots of x(Ca, r) calculated analytically from the scaling function, Eqs.(18) and (26) with parameter values from Table1.The fitting is quite consistent and accurate, both qualitatively and quantitatively provide an indicative direct comparison between plots of the source data (from DeProF model simulations) with plots of the scaling function derived in the previous section.The diagrams are paired.Diagrams of values predicted from DeProF model simulations are aligned on columns (a) and (c), Fig. 4. Plots of the analytical expression that fits the corresponding source data are aligned on columns (b) and (d).The fitting is quite consistent and accurate.

Fig. 7
Fig. 7 Energy maps, f EU (Ca, r) , based on Eqs.(30), for fluids with various n/w viscosity ratio values, .The maximum energy efficiency, f * EU (Ca, r * ) , as well as its horizontal projection, r * (Ca) , are depicted by the red, thick curves.The curtain of drop lines maps the locus of critical flow conditions, r * (Ca) ; it is suspended from a horizontal plane at a height equal to the nominal value of the maximum attainable energy efficiency for each two-fluid system (ceiling of efficiency), f * EU∞ .Thinner red curves are projections of f * EU (Ca, r * ) on the vertical planes.Ball markers are used to pinpoint order of magnitude changes in Ca 3024 on the top diagram, and f * EU∞ = 1∕(1 + 1.5) 2 = 0.160 on the bottom diagram.These limits are asymptotically reached at the corresponding nominal asymptotic values, r * ∞ , of the critical flow rate ratio at very high-capillary number (fully viscous flow, Ca → ∞ ), given by log r * ∞ = log � 1∕ √ � , i.e., log r * ∞ = −0.5 × log 0.67 = 0.0870 and log r * ∞ = −0.5 × log 1.5 = −0.0880, for the top and bottom diagram, respectively.

Fig. 8
Fig. 8 Universal, energy efficiency maps describing steady-state fully developed two-phase flow in pore networks for two values of n/w viscosity ratio, = 0.67 (top) and = 1.5 (bottom).The insert is the first impression of a universal diagram (adapted from Valavanides 2018b) C 3 y 3 + x − x 0 y 2 − C 0 = 0 (39) D = 4C 0 x − x 0 reduced cubic equation has a single root, given by Straightforward calculations lead to the following general expression,