© Pål Ø. Andersen, published by IFP Energies nouvelles, 2021

Nomenclature

Roman

C_s, C_i: Saturation and relative permeability slopes in the intercept method, −

F_i : Injected phase fraction, −

f_w : Water fractional flow function, −

J : Scaled capillary pressure, −

J₁, J₂, J₃: J-function coefficients, −

k_ri : Phase relative permeability, −

K : Absolute permeability, m²

k₁, k₂: J-function parameters, −

L : Core length, m

n : Exponent for approximating end effect saturation profile, −

n₁, n₂: J-function exponents, −

n_i : Phase Corey exponent, −

n_i1, n_i2: Phase Corey exponent end values, −

N₀ : capillary number (viscous to capillary forces), −

p_i : Phase pressure, Pa

P_c : Capillary pressure, Pa

s_i : Phase saturation, −

$s_w^{\mathrm{eq}}$ : Water saturation at which capillary pressure is zero, −

$s_w^r$ : Water saturation at a given fraction corrected for end effects, −

S_i : Normalized phase saturation, −

S₁ : Normalized water saturation at inlet (Y = 1), −

S_eq : Normalized water saturation at which capillary pressure is zero, −

S_r : Reference scaled saturation (obtained if no end effects present), −

$\overline{S}$ : Normalized water saturation averaged over the core, −

u_i : Darcy phase velocity, m/s

v_i : Interstitial velocity, m/s

Y : Scaled distance from outlet, −

Greek

ϕ : Porosity, −

μ_i : Phase viscosity, Pa s

σ_ow : Interfacial tension, N/m

λ_i : Phase mobility, 1/(Pa s)

Δp_i: Phase pressure drop, Pa

Δp_i,ref: Pressure drop without end effects, Pa

Indices

c : Condensate

eq: Zero capillary pressure condition

g : Gas

i : Phase index

o : Oil

r : Reference (no end effects)

T : Total

w : Water

1 Introduction

Relative permeabilities are an essential input to simulation when modeling multiphase flow processes, including petroleum recovery, carbon storage, and liquid invasion in hydraulic fracturing (Jeong et al., 2021; Juanes et al., 2006). Relative permeabilities describe the reduction of mobility for a flowing phase in presence of other phases and are modeled as function of saturation, although the saturation path, wettability, hysteresis, stress, and temperature and viscous coupling can affect them (Andersen et al., 2020a; Anderson, 1987; Bourbiaux and Kalaydjian, 1990). Their measurement should hence be performed under conditions close to the expected process taking place in the reservoir (Sidiq et al., 2017).

Relative permeabilities are commonly measured experimentally using the steady state technique, unsteady state technique or centrifuge. The steady state technique will be the focus of this work: It consists of injecting two fluids in different flow fractions (which can include fractions corresponding to single phase injection) and measuring the average saturation and the pressure drop over the core at steady state, when production rates and pressure readings have stabilized. If sufficient time has passed to reach steady state and the saturations are uniform, Darcy’s law can be applied to calculate relative permeabilities directly (Richardson et al., 1952). The unsteady state method considers injection of one phase to displace another where use is made of the transient production and pressure data before steady state is reached to calculate relative permeability (Johnson et al., 1959). With centrifuge, rotation at high speed allows calculating relative permeability for the less mobile phase from transient production data (Hagoort, 1980). All these methods are complicated by the presence of capillary pressure. This will be discussed only for the steady state method. Relative permeability can also be predicted from theoretical bundle of tubes models or digital rock models describing flow in representative porous structures (Nguyen et al., 2006; Valavanides, 2018).

In an open space there is no confinement or curvature of the fluid–fluid interface, setting capillary pressure equal zero at the producing end face of a core or an open fracture in a reservoir (Leverett, 1941). Pressure continuity forces capillary pressure profiles to converge to zero at the core outlet. Since there is a unique relation between capillary pressure and saturation (when saturation is changed monotonously) the steady state saturation profile will converge to the saturation giving zero capillary pressure. This results in nonuniform saturation profiles with most deviation at the outlet (Richardson et al., 1952). The capillary end effects can impact the calculation of relative permeabilities (Osoba et al., 1951). End effects cause steady state average saturation and pressure drop across the core to differ from what their values would be without end effects. Increasing the ratio of advective to capillary forces suppresses the end effect and makes the saturations more uniform and in accordance with the injected flow fraction. End effects can result in and explain apparently rate-dependent relative permeabilities and flooding behavior, as observed experimentally (Alizadeh et al., 2007; Henderson et al., 1998; Jeong et al., 2021; Odeh and Dotson, 1985; Osoba et al., 1951; Rapoport and Leas, 1953). Chen and Wood (2001) found rate-insensitive relative permeabilities when their in-situ measurements indicated negligible fluid accumulation at the outlet. Zou et al. (2020) used in-situ imaging of a core plug steady state saturation profile to calculate relative permeability, but required an independent measurement of the capillary pressure function. This was achieved both via experimental and pore scale imaging based techniques. Rapoport and Leas (1953) showed that by increasing the core length, injected fluid viscosity and injection rate the end effects were less significant. However, reaching sufficiently high rates that end effects are negligible is not always practical due to limitations on core integrity, flow rate capacity, pressure reading, flow regime, and unrepresentative mobilization of residual droplets.

By varying the flow rate at a given flow fraction it is possible to assess the importance of the end effects and correct the measurements to representative saturations and pressure drops. Gupta and Maloney (2016) found that plotting pressure drop against rate was linear and gave a constant pressure drop at zero rate associated with end effects. Removing this extra pressure drop at each rate gave corrected consistent relative permeabilities. Further, they linked the extent of the end effect profile (how much of the core it covered) to how large the added pressure drop was relative to the corrected pressure drop. Plotting average saturation against this ratio could be extrapolated linearly to when the end effect had zero extent. The saturation and relative permeabilities obtained from these linear extrapolations are corrected for end effects and based on where the lines intercept. They hence called this the intercept method. Their main assumption was that the end effect saturation profile was limited to within the core. No relation was made to physical system properties, flow conditions or saturation functions. Their assumptions were mainly justified by numerical and experimental examples, rather than theory. A review of the intercept method was given by Reed and Maas (2018).

Andersen et al. (2017a) derived explicit analytical solutions for water flooding with end effects assuming specific saturation function correlations. They derived the intercept method theoretically, expressed using that average saturation and pressure drop divided by rate behaved linearly with inverse rate. They found how the saturation and pressure line slopes were related to input parameters of saturation functions. A capillary number was found determining the linear behavior with a critical value of 1 determining when the linear behavior was valid. They also could predict behavior when the end effects were outside the core. The model was validated experimentally by Andersen et al. (2020b) to determine both water relative permeability and capillary pressure. Huang and Honarpour (1998) also considered the waterflooding case, but under different saturation function assumptions. Their solutions were implicit, but could be used to correct end point relative permeability. Andersen and Zhou (2020) considered co-injection tests under the constraint of linear saturation functions and derived a capillary number incorporating all system parameters including the saturation functions. Virnovsky et al. (1995) demonstrated how sensitivity in average saturation and phase pressure drops to injection rate were analytically related to saturation functions. Andersen et al. (2020b) and Santos et al. (2021) history matched relative permeability and capillary pressure from multi-rate tests. It has also been demonstrated that saturation functions can be determined from co-current spontaneous imbibition experiments by varying the viscosity ratio systematically (Andersen, 2021; Andersen et al., 2019).

In this work we present mathematical relations describing end effects during steady state tests from fundamental assumptions and derive general conclusions regarding saturation profiles, parameters that affect saturation profile shape or just compress the profile, the intercept method, scaling, and dimensionless numbers.

Calculating effective relative permeabilities with Darcy’s law from data with end effects will be inaccurate, but combining such data from different total rates can allow accurate prediction. We derive a more intuitive and more general intercept method compared to the work of Gupta and Maloney (2016) based on established assumptions in core scale multiphase flow simulation. We show, for given fluids and injected flow fraction but varied total rate, that plotting steady state average saturation and inverse effective relative permeability against inverse total rate yields straight lines (at high total rates). The lines result in correct relative permeability and saturation points for that flow fraction when extrapolated to zero inverse total rate (infinite total rate). It is shown how the slopes of these lines are related to the saturation functions, allowing us to scale all data universally to the same straight lines. Two dimensionless numbers are obtained that directly express when end effects become important and how much in terms of how much (a) the average saturation is changed and (b) the effective relative permeabilities are reduced compared to values unaffected by end effects. A third dimensionless number is derived with a critical value of ~1 as a universal criterion for when the linear trends with inverse total rate (the intercept method) begin. The critical value is directly reflecting that the end effect profile has reached the inlet end of the core. When the profile extends beyond the inlet, the trends become nonlinear. All the dimensionless numbers are divided into a part depending on saturation functions, injected flow fraction and viscosity ratio and a second part containing constant known fluid and rock parameters such as core length, porosity, interfacial tension, etc. Methodologies for deriving relative permeability and capillary pressure consistently, even based on combining data from tests with different fluid and core properties are presented and demonstrated on two datasets from the literature.

2. Theory

2.1 Transport equations

The mathematical description of 1D incompressible and immiscible flow of oil (o) and water (w) in a porous homogeneous medium under negligible influence of gravity is given by mass balance and Darcy’s law, respectively:$$\varphi \frac{\mathrm{\partial }{s}_i}{\mathrm{\partial }t}=-\frac{\mathrm{\partial }{u}_i}{\mathrm{\partial }x}, \hspace{1em}\left(i=o,w\right),$$(1)

$$u_i=-K{\lambda }_i\frac{\mathrm{\partial }{p}_i}{\mathrm{\partial }x}, \hspace{1em}{\lambda }_i=\frac{k_{\mathrm{ri}}}{{\mu }_i}, \hspace{1em}\left(i=o,w\right),$$(2)

ϕ is porosity, s_i saturation of phase i = o, w, u_i Darcy velocity, K absolute permeability, λ_i mobility, k_ri relative permeability, μ_i viscosity, and p_i pressure. The saturations are dependent due to volume conservation, and the pressures are related by the capillary pressure function:$$s_w+s_o=1, \hspace{1em}{p}_o-p_w=P_c\left(s_w\right).$$(3)

The total Darcy velocity u_T is defined as:$$u_T=u_o+u_w=-K{\lambda }_T\frac{\mathrm{\partial }{p}_w}{\mathrm{\partial }x}-K{\lambda }_o\frac{\mathrm{\partial }{P}_c}{\mathrm{\partial }x}.$$(4)

Also, λ_T is the total mobility, given by:$${\lambda }_T={\lambda }_o+{\lambda }_w.$$(5)

It follows from adding the transport equations in (1) that:$$\frac{\mathrm{\partial }{u}_T}{\mathrm{\partial }x}=0.$$(6)

The water phase equation can then be expressed with variables u_T, s_w as:$$\varphi \frac{\mathrm{\partial }{s}_w}{\mathrm{\partial }t}=-\frac{\mathrm{\partial }{u}_w}{\mathrm{\partial }x},$$(7)

$$u_w=u_T{f}_w+K{\lambda }_o{f}_w\frac{\mathrm{\partial }{P}_c}{\mathrm{\partial }x},$$(8)

where f_w is the fractional flow function defined by:$$f_w=\frac{{\lambda }_w}{{\lambda }_w+ {\lambda }_o}.$$(9)

2.2 Boundary and initial conditions

Water and oil are injected simultaneously at the inlet x = 0 with a water flow fraction F (the water fraction of the total injected flux) and a total Darcy flux u_T (Fig. 1):$$F=\frac{u_w\left(x=0\right)}{u_w\left(x=0\right)+u_o\left(x=0\right)}=\frac{u_w\left(x=0\right)}{u_T}.$$(10)

The injected water flux given u_T and F is then:$$u_w\left(x=0\right)=u_{T}F.$$(11)

The water flux (and that of oil) is composed of both an advective and capillary component, see (8). Hence, F does not correspond to f_w unless the capillary pressure gradient can be ignored. From (8) we write this boundary condition as:$$u_w\left(x=0\right)={\left[u_T{f}_w+K{\lambda }_o{f}_w\frac{\mathrm{\partial }{P}_c}{\partial x}\right]}_{x=0}=u_{T}F.$$(12)

The outlet boundary condition is described by a zero capillary pressure (Leverett, 1941), which corresponds to a fixed outlet water saturation:$$P_c\left(x=L\right)=0,\hspace{1em} s_w\left(x=L\right)=s_w^{\mathrm{eq}},$$(13)

where by definition P_c$\left(s_w^{\mathrm{eq}}\right)$ = 0.

2.3 Steady state

At steady state we have no changes with time in the system, i.e.:$$\frac{\mathrm{\partial }{s}_i}{\mathrm{\partial }t}=0, \hspace{1em}\frac{\mathrm{\partial }{p}_i}{\mathrm{\partial }t}=0, \hspace{1em}\left(i=o,w\right).$$(14)

The phases are non-uniformly distributed due to the balance between advective and capillary forces. Given that time is not influential at steady state; in the following, water saturation and water pressure will be taken as functions of spatial coordinate alone: s_w = s_w(x) and p_w = p_w(x). Equation (7) can then be written as:$$\frac{\mathrm{d}{u}_w}{\mathrm{d}x}=0=\frac{\mathrm{d}}{\mathrm{d}x}\left[u_T{f}_w+K\left({\lambda }_o{f}_w\frac{\mathrm{d}{P}_c}{\mathrm{d}x}\right)\right].$$(15)

At steady state the fluxes are uniform, i.e. the same amount of water and oil passes through every cross section, however the saturations and velocities can differ. Setting the water flux uniformly equal to that at the inlet, see (12), gives:$$u_w=u_{T}F=u_T{f}_w+K\left({\lambda }_o{f}_w\frac{\mathrm{d}{P}_c}{\mathrm{d}x}\right).$$(16)

Using that $\frac{\mathrm{d}{P}_c}{\mathrm{d}x}=\frac{\mathrm{d}{P}_c}{\mathrm{d}{s}_w}\frac{\mathrm{d}{s}_w}{\mathrm{d}x}$, we can solve (16) with respect to the saturation gradient:$$\frac{\mathrm{d}{s}_w}{\mathrm{d}x}=\frac{u_T\left(F-f_w\right)}{K{\lambda }_o{f}_w\frac{\mathrm{d}{P}_c}{\mathrm{d}{s}_w}}.$$(17)

The water saturation gradient is thus dependent on the two phase mobilities, the capillary pressure curve, the injected water flow fraction F and the injection flux u_T. We can further introduce the interstitial total velocity, v_T, and dimensionless Leverett J-function (Dullien, 2012):$$u_T=\varphi v_T, \hspace{1em}{P}_c={\sigma }_{\mathrm{ow}}\sqrt{\frac{\varphi }{K}}J\left(s_w\right),$$(18)

which results in:$$\frac{\mathrm{d}{s}_w}{\mathrm{d}x}=\frac{v_T\sqrt{\frac{\varphi }{K}}}{{\sigma }_{\mathrm{ow}}}\frac{\left(F-f_w\right)}{f_w{\lambda }_o \frac{\mathrm{d}J}{\mathrm{d}{s}_w}}.$$(19)

The above equation can be integrated to find the saturation distribution starting from s_w(x = L) = $s_w^{\mathrm{eq}}$. The saturation gradient will be nonzero until a saturation $s_w^r$ is reached such that:$$f_w\left(s_w^r\right)=F.$$(20)

After which the saturation remains stable at $s_w^r$. This is the state corresponding to negligible end effects. $s_w^r$ is found by solving (20). The pressure gradients of oil and water at steady state follow from (2) combined with (16):$$\frac{\mathrm{d}{p}_w}{\mathrm{d}x}=-\frac{u_{T}F}{K{\lambda }_w},\hspace{1em}\frac{\mathrm{d}{p}_o}{\mathrm{d}x}=-\frac{u_T\left(1-F\right)}{K{\lambda }_o}.$$(21)

The above corresponds to Darcy’s law, where the water and oil fluxes are constant equal to u_TF and u_T(1 − F) and the mobilities vary along the core according to the steady state saturation distribution found from (19).

2.4 Scaled saturation profile

2.4.1 Derivation

Assume a domain where every saturation has a unique position. Equation (19) can then be solved by separation into a space coordinate integral and a saturation integral:$$\frac{v_T\sqrt{\frac{\varphi }{K}}}{{\sigma }_{\mathrm{ow}}}{\int }_{x^{\text{'}}=L}^x{\mathrm{d}x}^{\mathrm{\prime }}={\int }_{s_w^{\mathrm{eq}}}^{s_w}\frac{\left[f_w{\lambda }_o\frac{\mathrm{d}J}{\mathrm{d}{s}_w}\right]\left(s_w^{\mathrm{\prime }}\right)}{F-f_w\left(s_w^{\mathrm{\prime }}\right)}{\mathrm{d}s}_w^{\mathrm{\prime }}={\int }_{S_{\mathrm{eq}}}^S\frac{\left[f_w{\lambda }_o\frac{\mathrm{d}J}{\mathrm{d}S}\right]\left({\mathrm{S}}^{\mathrm{\prime }}\right)}{F-f_w\left({\mathrm{S}}^{\mathrm{\prime }}\right)}{\mathrm{d}S}^{\mathrm{\prime }}.$$(22)

Although the former integral is trivial, the latter in most cases requires numerical methods. Note that the latter saturation integral above has been expressed using normalized saturation S which provides the following relations:$$S=\frac{s_w-s_{\mathrm{wr}}}{\Delta s_w}, \hspace{1em}\Delta s_w=1-s_{\mathrm{or}}-s_{\mathrm{wr}},$$(23)

$$\mathrm{d}{s}_w=\frac{\mathrm{d}{s}_w}{\mathrm{d}S}\mathrm{d}S=\Delta s_w\mathrm{d}S, \hspace{1em}{S}_{\mathrm{eq}}=\frac{s_w^{\mathrm{eq}}-s_{\mathrm{wr}}}{\Delta s_w},$$(24)

$$\frac{\mathrm{d}J}{\mathrm{d}{s}_w}=\frac{\mathrm{d}J}{\mathrm{d}S}\frac{\mathrm{d}S}{\mathrm{d}{s}_w}.$$(25)

The parameters s_wr, s_or denote the critical saturations of water and oil, respectively, where their respective relative permeability is zero. Δs_w denotes the magnitude of the mobile saturation interval. S_eq is the normalized saturation where capillary pressure is zero.

Note that the terms f_w and f_wλ_o can both be written as functions of viscosity ratio, and that the latter term is inversely proportional to the geometric mean of viscosities:$$f_w=\frac{k_{\mathrm{rw}}}{k_{\mathrm{rw}}+\left(\frac{{\mu }_w}{{\mu }_o}\right)k_{\mathrm{ro}}},\hspace{1em}{f}_w{\lambda }_o=\frac{1}{{\left({\mu }_o{\mu }_w\right)}^{0.5}}\frac{k_{\mathrm{rw}}{k}_{\mathrm{ro}}{\left(\frac{{\mu }_w}{{\mu }_o}\right)}^{0.5}}{k_{\mathrm{rw}}+\left(\frac{{\mu }_w}{{\mu }_o}\right)k_{\mathrm{ro}}}.$$(26)

From this it is convenient to introduce the notations:$$Y=\frac{y}{L}=\frac{L-x}{L},\hspace{1em}{N}_0=\frac{v_{T}L\sqrt{\frac{\varphi }{K}}{\mu }_m}{{\sigma }_{\mathrm{ow}}}, \hspace{1em}{\mu }_m={\left({\mu }_o{\mu }_w\right)}^{0.5}.$$(27)

N₀ is a dimensionless capillary number (ratio of viscous to capillary forces) containing static or single phase flow parameters. This number is fixed if the same fluids, core and total rate are considered. N₀ > 0 since v_T > 0 (flow in x-direction towards outlet). This leads to the solution form of interest:$$Y\left(S\right)=-\frac{1}{N_0}{\int }_{S_{\mathrm{eq}}}^S\frac{{\mu }_m\left[f_w{\lambda }_o\frac{\mathrm{d}J}{\mathrm{d}S}\right]\left({\mathrm{S}}^{\mathrm{\prime }}\right)}{F-f_w\left({\mathrm{S}}^{\mathrm{\prime }}\right)}{\mathrm{d}S}^{\mathrm{\prime }},$$(28)

which is valid for all S between S_eq and S_r. S_r is the normalized saturation obtained at the flow fraction F without end effects:$$S_r=\frac{s_w^{\mathrm{r}}-s_{\mathrm{wr}}}{\Delta s_w}.$$(29)

2.4.2 Shape characteristics

Consider two saturations S₁, S₂ with positions given by (28). Their relative position is given by:$$\frac{Y\left(S_2\right)}{Y\left(S_1\right)}=\frac{{\int }_{S_{\mathrm{eq}}}^{S_2}\frac{\left[f_w{\lambda }_o\frac{\mathrm{d}J}{\mathrm{d}S}\right]\left({\mathrm{S}}^{\mathrm{\prime }}\right)}{F-f_w\left({\mathrm{S}}^{\mathrm{\prime }}\right)}{\mathrm{d}S}^{\mathrm{\prime }}}{{\int }_{S_{\mathrm{eq}}}^{S_1}\frac{\left[f_w{\lambda }_o\frac{\mathrm{d}J}{\mathrm{d}S}\right]\left({\mathrm{S}}^{\mathrm{\prime }}\right)}{F-f_w\left({\mathrm{S}}^{\mathrm{\prime }}\right)}{\mathrm{d}S}^{\mathrm{\prime }}}.$$(30)

This indicates that only the saturation functions, viscosity ratio and injected fraction determine the shape. The parameters in N₀ will affect how compressed the profile is, but not its shape.

2.4.3 Length and range of end effect region

Based on the previous notation (27) the saturation gradient can be written from (19) as:$$\frac{\mathrm{d}S}{\mathrm{d}Y}=-N_0\frac{F-f_w}{{\mu }_m{f}_w{\lambda }_o\frac{\mathrm{d}J}{\mathrm{d}S}}=\alpha \left(S\right)\left(F-f_w\right).$$(31)

Assume first that two phases are injected simultaneously (0 < F < 1) which means S_r ∉ {0, 1}. Further, consider saturations far from the outlet Y ≫ 0. The coefficient $\alpha \left(S\right)=-\frac{N_0}{{\mu }_m{f}_w{\lambda }_o\frac{\mathrm{d}J}{\mathrm{d}S}}>0$ will then obey 0 < α₀ < α(S) < α₁ for some finite limits α₀, α₁. As distance from outlet increases the saturation will approach S_r according to the following Taylor approximation:$$\mathrm{d}Y=\frac{1}{\alpha \left(S_r\right)}\frac{\mathrm{d}S}{\left(F-f_w\right)}\approx \frac{1}{\alpha \left(S_r\right)\frac{\mathrm{d}{f}_w}{\mathrm{d}S}{|}_{S_r}}\frac{\mathrm{d}\mathrm{S}}{\left(S_r-S\right)}.$$(32)

Y will increase when S → S_r. We have considered a saturation region sufficiently close to S_r that α and $\frac{\mathrm{d}{f}_w}{\mathrm{d}S}$ can be considered constant as if evaluated at S_r. Assume a saturation S^* in that region with a finite difference from S_r. We integrate (32) and find that the spatial distance between the two saturations is infinite:$$Y\left(S_r\right)-Y\left(S^{\mathrm{*}}\right)=-\frac{1}{\alpha \left(S_r\right)\frac{\mathrm{d}{f}_w}{\mathrm{d}S}{|}_{S_r}}\underset{\mathrm{S}\to S_r}{\lim }\ln \left(\frac{S_r-S}{S_r-S^{*}}\right)=\mathrm{\infty }.$$(33)

S_r is only reached at infinite distance, in other words, it takes an infinite distance for end effects to vanish. This was exemplified by Andersen and Zhou (2020) for the special case of linear saturation functions, but is now proved in general.

On the other hand, assume single phase injection as given by F ∈ {0, 1}, which also implies S_r ∈ {0, 1}. The term f_wλ_o$\frac{\mathrm{d}J}{\mathrm{d}S}$ is the saturation dependent part of what is referred to as the capillary diffusion coefficient. It has parabolic shape and equals zero (only) at the points S = {0, 1}. Assume scaled saturations S approaching S_r and evaluate the saturation gradient:$$\frac{\mathrm{d}S}{\mathrm{d}Y}{\mathrm{}|}_{S_r}=-N_0{\mu }_m^{-1}\underset{\mathrm{S}\to S_r}{\lim }\frac{\left(F-f_w\right)}{f_w{\lambda }_o\frac{\mathrm{d}J}{\mathrm{d}S}}=N_0{\mu }_m^{-1}\underset{S\to S_r}{\lim }\frac{\frac{\mathrm{d}{f}_w}{\mathrm{d}S}}{\frac{\mathrm{d}}{\mathrm{d}S}\left(f_w{\lambda }_o\frac{\mathrm{d}J}{\mathrm{d}S}\right)} .$$(34)

As both nominator and denominator in the first limit approach zero, we have applied L’Hopital’s rule. There exist parameter choices where $\frac{\mathrm{d}{f}_w}{\mathrm{d}S}$ could be zero or nonzero at the end points; using Corey relative permeabilities with exponents greater than 1 or equal to 1, respectively, will give such behavior Andersen et al. (2017a, 2020b) showed that the use of Corey-type functions for both relative permeability and capillary pressure resulted in a finite length end effect. Huang and Honarpour (1998) on the other hand used Corey-Burdine equations and obtained an infinite length.

The above derivations and examples show that the end effect will always have infinite length when two fluids are injected simultaneously. During single phase injection, the end effect will have either infinite or finite length depending on the saturation function correlations. We will however show that regardless of injection conditions and saturation function correlations we can define a practical length which produces the same results as if there was a finite length, and thus the intercept method can be derived. The intercept method was originally derived based entirely on the assumption that the end effect did not exceed the length of the core (Gupta and Maloney, 2016). It was found from several numerical examples that although the position of the saturation S_r extends to infinity, key saturation integrals converged as the integral limit approached S_r.

We have demonstrated the validity of the continuous saturation profile (28) from S_eq to S_r where a unique relation exists between Y and S. However, if S_r is finally reached at a specific Y(S_r), as can happen for the single phase injection case, all Y > Y(S_r) will have a constant saturation S_r which is necessary to account for as it impacts calculations of average saturation and pressure drop.

2.5 Dimensionless numbers based on cumulative end effect profile

Consider for simplicity that the porous medium extends infinitely beyond the core length Y = 1 such that saturations between S_eq and S_r have specified positions Y(S). The area between the straight line S = S_r and the graph Y(S) represents the cumulative amount of phase trapped by end effects as measured in displaceable pore volumes, n_dpv:$$n_{\mathrm{dpv}}={\int }_{S=S_{\mathrm{eq}}}^{S_r}Y\left|\mathrm{d}S\right|=-\frac{1}{N_0}{\int }_{S=S_{\mathrm{eq}}}^{S_r}{\int }_{S_{\mathrm{eq}}}^S\frac{{\mu }_m\left[f_w{\lambda }_o\frac{\mathrm{d}J}{\mathrm{d}S}\right]\left({\mathrm{S}}^{\mathrm{\prime }}\right)}{F-f_w\left({\mathrm{S}}^{\mathrm{\prime }}\right)}{\mathrm{d}S}^{\mathrm{\prime }}\left|\mathrm{d}S\right|.$$(35)

The absolute sign on dS in the outer integral is used to produce a positive value regardless of whether S_eq is larger or less than S_r. As an example, if S_r = 0, S_eq = 1 and Y(S) = 1 (all the saturations 0 < S < 1 are positioned at the inlet end of the core), the area is 1 meaning n_dpv = 1. Whether the end effects have impact on the system will depend on how far the saturation profile is deviated from S_r throughout the core. n_dpv unfortunately cannot distinguish whether the accumulated phase is located within the core or not, but is a strong indicator of impact if the end effects are known to be within the core. The average length Y_av of the end effect profile (in core lengths) is:$$Y_{\mathrm{av}}=\frac{1}{S_r-S_{\mathrm{eq}}}{\int }_{S=S_{\mathrm{eq}}}^{S_r}Y\mathrm{d}S=-\frac{1}{S_r-S_{\mathrm{eq}}}\frac{1}{N_0}{\int }_{S=S_{\mathrm{eq}}}^{S_r}{\int }_{S_{\mathrm{eq}}}^S\frac{{\mu }_m\left[f_w{\lambda }_o\frac{\mathrm{d}J}{\mathrm{d}S}\right]\left({\mathrm{S}}^{\mathrm{\prime }}\right)}{F-f_w\left({\mathrm{S}}^{\mathrm{\prime }}\right)}{\mathrm{d}S}^{\mathrm{\prime }}\mathrm{d}S.$$(36)

If Y_av = 0 there are no end effects while if Y_av = 1 the end effects are sufficiently strong to have a profile with average distance equal the length of the core. The saturation profile is nonuniform. Especially, the profile is locked at Y = 0 at S_eq so the saturations closer to S_r have distances greater than Y_av. Assume the profile Y(S) is approximated by an n’th order polynomial $\tilde{Y}$ of S, such that Y_av is preserved and $\tilde{Y}$(S_eq) = 0:$$\tilde{Y}\left(S\right)=Y_{\mathrm{av}}\left(n+1\right){\left(\frac{S-S_{\mathrm{eq}}}{S_r-S_{\mathrm{eq}}}\right)}^n.$$(37)

For n = 1 the profile is linear while larger n give the profile more curvature and results in the end effect region to terminate at position Y_av(n + 1):$$Y_{\mathrm{cee}}=\tilde{Y}\left(S_r\right)=Y_{\mathrm{av}}\left(n+1\right).$$(38)

The parameter n can be selected to fit typical profiles or given an assumed value, but will be assumed fixed.

n_dpv and Y_cee are dimensionless capillary numbers expressing in different forms the ratio of capillary to viscous forces. They are derived from physical considerations and uniquely combine any set of system parameters, including the saturation functions. Large values ≫1 indicate strong end effects, low values ≪1 indicate negligible end effects and values of magnitude ≈1 are expected to give a transition in behavior. Especially, when Y_cee ≈ 1 the end effect profile should exactly reach the inlet. The two numbers are related by:$$Y_{\mathrm{cee}}=\frac{n+1}{\left|S_{\mathrm{eq}}-S_r\right|}{n}_{\mathrm{dpv}}\gg n_{\mathrm{dpv}}.$$(39)

Whether the end effects are limited to within the core or not is essential to the derivation of the intercept method (Andersen and Zhou, 2020; Andersen et al., 2017a, 2020b; Gupta and Maloney, 2016), suggesting Y_cee = 1 to be a criterion for when the method is valid.

2.6 Average saturation

2.6.1 General definition

In the following it will be assumed that the saturation profile extends to infinity, i.e. either co-injection is considered or we have single phase injection with a proper combination of saturation functions. The following spatial integrals are independent of this assumption, but their conversion to saturation integrals are not.

The core average saturation follows from integrating the saturation along the core. The integral can be converted into a saturation integral evaluated from the scaled outlet saturation S^eq to the scaled saturation S₁ at the inlet Y = 1:$$\overline{S}={\int }_{Y=0}^{1}S\left(Y\right)\mathrm{d}Y={\int }_{S=S_{\mathrm{eq}}}^{S_1}S\frac{\mathrm{d}Y}{\mathrm{d}S}\mathrm{d}S=-\frac{1}{N_0}{\int }_{S=S_{\mathrm{eq}}}^{S_1}S\frac{{\mu }_m{f}_w{\lambda }_o\frac{\mathrm{d}J}{\mathrm{d}S}}{F-f_w}\mathrm{d}S.$$(40)

In the above we have applied:$$\frac{\mathrm{d}Y}{\mathrm{d}S}=-\frac{1}{N_0}\frac{{\mu }_m{f}_w{\lambda }_o\frac{\mathrm{d}J}{\mathrm{d}S}}{F-f_w},$$(41)

which follows from (31). The inlet saturation S₁ is unknown and depends on the extent of the saturation profile. It is found by solving the equation Y(S₁) = 1, equivalently:$$1=-\frac{1}{N_0}{\int }_{S_{\mathrm{eq}}}^{S_1}\frac{{\mu }_m\left[f_w{\lambda }_o\frac{\mathrm{d}J}{\mathrm{d}S}\right]\left(S^{\mathrm{\prime }}\right)}{F-f_w}{\mathrm{d}S}^{\mathrm{\prime }}.$$(42)

As seen, S₁ depends on N₀, S_eq and the saturation functions. Especially, the saturation at this or any other specific location will depend on both the general profile shape and how compressed it is (the magnitude of N₀).

2.6.2 Intercept method

Assume now the saturation profile is approximated by introducing a saturation S^* such that:$$S^{\mathrm{*}}=S_r\left(1-\epsilon \right)+S_{\mathrm{eq}}\epsilon , \hspace{1em}0<\epsilon <1,$$(43)

for some small ε where:$$Y\left(S\right)=-\frac{1}{N_0}{\int }_{S_{\mathrm{eq}}}^S\frac{{\mu }_m\left[f_w{\lambda }_o\frac{\mathrm{d}J}{\mathrm{d}S}\right]\left({\mathrm{S}}^{\mathrm{\prime }}\right)}{F-f_w\left({\mathrm{S}}^{\mathrm{\prime }}\right)}{\mathrm{d}S}^{\mathrm{\prime }},\hspace{1em}\left(Y<Y\left(S^{\mathrm{*}}\right)\right),$$(44)

$$S=S_r, \hspace{1em}\left(Y>Y\left(S^{\mathrm{*}}\right)\right).$$(45)

In other words, at a given fraction (1 − ε) of the saturation interval between S_eq and S_r, at a saturation S^* sufficiently close to S_r: the (infinite) end effect saturation profile is approximated to a finite end effect saturation profile stopping at Y(S^*). At greater distances there are assumed no end effects, S = S_r. The closer S^* is to S_r the better is the approximation. Further, assume that the length of the end effect is within the core: Y(S^*) < 1. This distance will be called Y^*. There is then a corresponding distance 1 − Y^* at the inlet without end effects. The average saturation in the core under such circumstances can be described by:$$\overline{S}={\overline{S}}^{*}{Y}^{*}+S_r\left(1-Y^{*}\right),$$(46)

$$Y^{\mathrm{*}}=-\frac{1}{N_0}{\int }_{S_{\mathrm{eq}}}^{S^{\mathrm{*}}}\frac{{\mu }_m\left[f_w{\lambda }_o\frac{\mathrm{d}J}{\mathrm{d}S}\right]\left({\mathrm{S}}^{\mathrm{\prime }}\right)}{F-f_w\left({\mathrm{S}}^{\mathrm{\prime }}\right)}{\mathrm{d}S}^{\mathrm{\prime }}.$$(47)

${\overline{S}}^{\mathrm{*}}$ is the average saturation in the end effect region (0 < Y < Y^*) and the position Y^* is known. By considering the saturation profile between S_eq and S^* we find the average saturation in that interval ${\overline{S}}^{\mathrm{*}}$:$${\overline{S}}^{\mathrm{*}}=\frac{1}{Y^{\mathrm{*}}}{\int }_{Y=0}^{Y^{\mathrm{*}}}S\left(Y\right)\mathrm{d}Y=\frac{1}{Y^{\mathrm{*}}}{\int }_{S=S_{\mathrm{eq}}}^{S^{\mathrm{*}}}S\frac{\mathrm{d}Y}{\mathrm{d}S}\mathrm{d}S=\frac{{\int }_{S_{\mathrm{eq}}}^{S^{\mathrm{*}}}S\frac{{\mu }_m\left[f_w{\lambda }_o\frac{\mathrm{d}J}{\mathrm{d}S}\right]\left(S\right)}{F-f_w}\mathrm{d}S}{{\int }_{S_{\mathrm{eq}}}^{S^{\mathrm{*}}}\frac{{\mu }_m\left[f_w{\lambda }_o\frac{\mathrm{d}J}{\mathrm{d}S}\right]\left({\mathrm{S}}^{\mathrm{\prime }}\right)}{F-f_w\left({\mathrm{S}}^{\mathrm{\prime }}\right)}{\mathrm{d}S}^{\mathrm{\prime }}}.$$(48)

As seen the average saturation ${\overline{S}}^{\mathrm{*}}$ in the end effect region is only a function of the saturation functions, flow fraction and viscosity ratio and not the parameters in N₀. The average saturation $\overline{S}$ in the core however varies linearly with Y^*, the fraction of the core covered by end effects. Gupta and Maloney (2016) assumed this and verified it by running numerical simulations, but did not make a formal proof.

From the definition of Y^* and ${\overline{S}}^{\mathrm{*}}$ in (47) and (48) we can write (46) as:$$\overline{S}=S_r+\frac{1}{N_0}{\int }_{S=S_{\mathrm{eq}}}^{S^{\mathrm{*}}}\left(S_r-S\right)\frac{{\mu }_m{f}_w{\lambda }_o\frac{\mathrm{d}J}{\mathrm{d}S}}{F-f_w}\mathrm{d}S,$$(49)

showing that the average saturation is linear with the inverse capillary number where the slope is a saturation integral and the intercept is the saturation without end effects.

The combination of the saturation integrals in (47) and (48) to one makes the final integral in (49) less sensitive to the choice of S^* and we can in fact let S^* → S_r. We then get a more correct (and less subjective) slope since at higher capillary numbers more of the profile will be within the core and should be accounted for.

The average (absolute) saturation can then be expressed as linear with the inverse capillary number 1/N₀ with a saturation term slope C_s and, by expanding the capillary number expression, proportional with inverse velocity:$${\overline{s}}_w=s_w^r+\frac{C_s}{N_0}=s_w^r+C_s\frac{{\sigma }_{\mathrm{ow}}}{L\sqrt{\frac{\varphi }{K}}{\mu }_m}\frac{1}{v_T},$$(50)

$$C_s=\Delta s_w{\int }_{S=S_{\mathrm{eq}}}^{S_r}\left(S_r-S\right)\frac{{\mu }_m{f}_w{\lambda }_o\frac{\mathrm{d}J}{\mathrm{d}S}}{F-f_w}\mathrm{d}S={\int }_{s_w=s_w^{\mathrm{eq}}}^{s_w^r}\left(s_w^r-s_w\right)\frac{{\mu }_m{f}_w{\lambda }_o\frac{\mathrm{d}J}{\mathrm{d}{s}_w} }{F-f_w}\mathrm{d}{s}_w.$$(51)

The intercept is the corrected saturation without end effects $s_w^r$.

2.7 Pressure analyses

2.7.1 Gradients

The pressure gradients of oil and water can be expressed in terms of the scaled distance from the outlet:$$\frac{\mathrm{d}{p}_w}{\mathrm{d}Y}=\frac{\mathrm{L\varphi }{v}_T{\mu }_{w}F}{K{k}_{\mathrm{rw}}\left(S\left(Y\right)\right)}>0,$$(52)

$$\frac{\mathrm{d}{p}_o}{\mathrm{d}Y}=\frac{\mathrm{L\varphi }{v}_T{\mu }_o\left(1-F\right)}{K{k}_{\mathrm{ro}}\left(S\left(Y\right)\right)}>0.$$(53)

2.7.2 Pressure drop

We define the pressure drop of a phase as the pressure at the inlet minus that at the outlet. We obtain these parameters by integration of the pressure gradients, expressed either as integrals over the positions of the core or saturation integrals, where we make use of knowing the saturation profiles:$$\Delta p_w=-\frac{1}{N_0}\frac{\mathrm{L\varphi }{v}_T{\mu }_m}{K}F{\int }_{S^{\mathrm{eq}}}^{S_1}\frac{1-f_w}{F-f_w}\frac{\mathrm{d}J}{\mathrm{d}S}\mathrm{d}S,$$(54)

$$\Delta p_o=-\frac{1}{N_0}\frac{\mathrm{L\varphi }{v}_T{\mu }_m}{K}\left(1-F\right){\int }_{S^{\mathrm{eq}}}^{S_1}\frac{f_w}{F-f_w}\frac{\mathrm{d}J}{\mathrm{d}S}\mathrm{d}S.$$(55)

Assume for comparison that there were no end effects such that the pressure gradients (52) and (53) are constant and evaluated at S_r:$${\left(\frac{\mathrm{d}{p}_w}{\mathrm{d}Y}\right)}_{\mathrm{ref}}=\frac{\mathrm{L\varphi }{v}_T{\mu }_{w}F}{K{k}_{\mathrm{rw}}\left(S_{\mathrm{r}}\right)}>0,$$(56)

$${\left(\frac{\mathrm{d}{p}_o}{\mathrm{d}Y}\right)}_{\mathrm{ref}}=\frac{\mathrm{L\varphi }{v}_T{\mu }_o\left(1-F\right)}{K{k}_{\mathrm{ro}}\left(S_r\right)}>0.$$(57)

Integrating from Y = 0 to 1, the corresponding pressure drops are:$$\Delta p_{w,\mathrm{ref}}=\frac{\mathrm{L\varphi }{v}_T{\mu }_{w}F}{K{k}_{\mathrm{rw}}\left(S_{\mathrm{r}}\right)}>0,$$(58)

$$\Delta p_{o,\mathrm{ref}}=\frac{\mathrm{L\varphi }{v}_T{\mu }_o\left(1-F\right)}{K{k}_{\mathrm{ro}}\left(S_r\right)}>0.$$(59)

From the definition of S_r we have that f_w(S_r) = F which can be expanded to:$$\frac{k_{\mathrm{rw}}\left(S_r\right)}{{\mu }_w{\lambda }_T\left(S_r\right)}=F,$$(60)

and equivalently expressed for the oil phase using 1 − f_w(S_r) = 1 − F, we have:$$\frac{k_{\mathrm{ro}}\left(S_r\right)}{{\mu }_o{\lambda }_T\left(S_r\right)}=1-F.$$(61)

The two equations are both related to the total mobility and can therefore be combined:$$\frac{{\mu }_{w}F}{k_{\mathrm{rw}}\left(S_{\mathrm{r}}\right)}=\frac{{\mu }_o\left(1-F\right)}{k_{\mathrm{ro}}\left(S_{\mathrm{r}}\right)}=\frac{1}{{\lambda }_T\left(S_{\mathrm{r}}\right)}.$$(62)

We thus see that the phase dependent terms in (56)–(59) are equal which implies that the phases have identical pressure gradients and pressure drop over the core in absence of end effects:$${\left(\frac{\mathrm{d}{p}_i}{\mathrm{d}Y}\right)}_{\mathrm{ref}}={\left(\frac{\mathrm{d}p}{\mathrm{d}Y}\right)}_{\mathrm{ref}}=\frac{\mathrm{L\varphi }{v}_T}{K{\lambda }_T\left(S_{\mathrm{r}}\right)},$$(63)

$$\Delta p_{i,\mathrm{ref}}=\Delta p_{\mathrm{ref}}=\frac{\mathrm{L\varphi }{v}_T}{K{\lambda }_T\left(S_{\mathrm{r}}\right)}.$$(64)

We then also obtain identical phase pressure gradients once considering positions adequately far from the end effect zone. The ratios of pressure drop with end effects to pressure drop without end effects are:$$\frac{\Delta p_w}{\Delta p_{\mathrm{ref}}}=-\frac{1}{N_0}\frac{k_{\mathrm{rw}}\left(S_r\right){\mu }_m}{{\mu }_w}{\int }_{S_{\mathrm{eq}}}^{S_1}\frac{1-f_w}{F-f_w}\frac{\mathrm{d}J}{\mathrm{d}S}\mathrm{d}S>0,$$(65)

$$\frac{\Delta p_o}{\Delta p_{\mathrm{ref}}}=-\frac{1}{N_0}\frac{k_{\mathrm{ro}}\left(S_r\right){\mu }_m}{{\mu }_o}{\int }_{S_{\mathrm{eq}}}^{S_1}\frac{f_w}{F-f_w}\frac{\mathrm{d}J}{\mathrm{d}S}\mathrm{d}S>0.$$(66)

It is seen that N₀ controls the impact of end effects on pressure drop for each phase, with the exception that the geometric viscosity in N₀ is replaced by the phase viscosity (since N₀ ∝ μ_m).

Of most interest from the pressure measurements are the data corresponding to the relative permeabilities without end effects k_ri(S_r). Assume that an “effective relative permeability” ${\tilde{k}}_{\mathrm{ri}}$ is calculated based on the measured pressure drop and injection conditions by direct application of Darcy’s law. On the other hand, the “true” relative permeability without end effects k_ri(S_r) would be obtained if the pressure drop was Δp_ref. This is related as follows:$$k_{\mathrm{ri}}\left(S_{\mathrm{r}}\right)=\frac{\mathrm{L\varphi }{v}_T{\mu }_i{F}_i}{K\Delta p_{\mathrm{ref}}}, \hspace{1em}{\tilde{k}}_{\mathrm{ri}}=\frac{\mathrm{L\varphi }{v}_T{\mu }_i{F}_i}{K\Delta p_i}.$$(67)

The ratio of pressure drops is then directly related to the ratio of relative permeability estimates:$$\frac{\Delta p_i}{\Delta p_{\mathrm{ref}}}=\frac{k_{\mathrm{ri}}\left(S_{\mathrm{r}}\right)}{{\tilde{k}}_{\mathrm{ri}}}.$$(68)

This implies the effective relative permeability approaches the true relative permeability when Δp_i approaches Δp_ref at high rates (small end effects):$${\tilde{k}}_{\mathrm{ri}}=k_{\mathrm{ri}}\left(S_{\mathrm{r}}\right)\frac{\Delta p_{\mathrm{ref}}}{\Delta p_i}.$$(69)

If the pressure drop is higher with end effects $\left(\frac{\Delta p_{\mathrm{ref}}}{\Delta p_i}<1\right)$ this leads to underestimation of the relative permeability $\left({\tilde{k}}_{\mathrm{ri}}<k_{\mathrm{ri}}\left(S_{\mathrm{r}}\right)\right)$, and vice versa.

2.7.3 The challenge of missing phase pressure data and its solution

Under normal circumstances we do not have access to the pressure drop of both phases, but measure only one phase pressure drop across the core. Virnovsky et al. (1995) state that during imbibition injection the pressures at the outlet are continuous, however only the nonwetting phase is continuous with the surroundings at the inlet. The experimentally measured pressure drop is therefore that of the nonwetting phase. Special inlet designs allowed pressure continuity and measurement of both fluid pressures (Virnovsky et al., 1998), but they are not common. An important question is then whether the available pressure drop still gives meaningful data to calculate the other phase’s relative permeability. Assume therefore that for a given phase the effective relative permeability is calculated using the pressure drop of the other phase. Equation (67) is then modified to:$$k_{\mathrm{ri}}\left(S_{\mathrm{r}}\right)=\frac{\mathrm{L\varphi }{v}_T{\mu }_i{F}_i}{K\Delta p_{\mathrm{ref}}}, \hspace{1em}{\tilde{k}}_{\mathrm{ri}}=\frac{\mathrm{L\varphi }{v}_T{\mu }_i{F}_i}{K\Delta p_j}, \hspace{1em}\left(j\ne i\in o,w\right),$$(70)

such that:$$\frac{k_{\mathrm{ri}}\left(S_{\mathrm{r}}\right)}{{\tilde{k}}_{\mathrm{ri}}}=\frac{\Delta p_j}{\Delta p_{\mathrm{ref}}}, \hspace{1em}\left(j\ne i\in o,w\right).$$(71)

Since Δp_j approaches Δp_ref when end effects are negligible the correct relative permeability is still obtained, although the effective estimates with end effects are different when using the pressure drop of the other phase.

2.7.4 Intercept method

Again, assume the saturation profile is divided into a region 0 < Y < Y^* with end effects and a region Y^* < Y < 1 without end effects where S = S_r, as described before, where S(Y^*) ≈ S_r. Particularly, the end effect region is within the core and there exists a region without end effects. By integrating the pressure gradients over each interval, the pressure drops of each phase are:$$\Delta p_w=\frac{\mathrm{L\varphi }{v}_T{\mu }_{w}F}{K}\left[-\frac{{\mu }_m}{N_0{\mu }_w}{\int }_{S=S_{\mathrm{eq}}}^{S^{\mathrm{*}}}\frac{1-f_w}{F-f_w}\frac{\mathrm{d}J}{\mathrm{d}S}\mathrm{d}S+\frac{1-Y^{\mathrm{*}}}{k_{\mathrm{rw}}\left(S_r\right)}\right],$$(72)

$$\Delta p_o=\frac{\mathrm{L\varphi }{v}_T{\mu }_o\left(1-F\right)}{K}\left[-\frac{{\mu }_m}{N_0{\mu }_o}{\int }_{S=S_{\mathrm{eq}}}^{S^{\mathrm{*}}}\frac{f_w}{F-f_w}\frac{\mathrm{d}J}{\mathrm{d}S}\mathrm{d}S+\frac{1-Y^{\mathrm{*}}}{k_{\mathrm{ro}}\left(S_r\right)}\right].$$(73)

Dividing the pressure drops Δp_i,ref with the reference pressure drops Δp_i (obtained with no end effects) and using the definition of Y^* from (47) we can express the result in the following form:$$\frac{\Delta p_w}{\Delta p_{w,\mathrm{ref}}}= 1+\frac{1}{N_0}\frac{{\mu }_m}{{\mu }_w}{\int }_{S=S_{\mathrm{eq}}}^{S^{\mathrm{*}}}\left[k_{\mathrm{rw}}-k_{\mathrm{rw}}\left(S_{\mathrm{r}}\right)\right]\frac{\left(1-f_w\right)\frac{\mathrm{d}J}{\mathrm{d}S}}{F-f_w} \mathrm{d}S,$$(74)

$$\frac{\Delta p_o}{\Delta p_{o,\mathrm{ref}}}= 1+\frac{1}{N_0}\frac{{\mu }_m}{{\mu }_o}{\int }_{S_{\mathrm{eq}}}^{S^{\mathrm{*}}}\left[k_{\mathrm{ro}}-k_{\mathrm{ro}}\left(S_r\right)\right]\frac{f_w}{F-f_w}\frac{\mathrm{d}J}{\mathrm{d}S}\mathrm{d}S.$$(75)

Again we note that the updated combined saturation integrals are less sensitive to S^* and we can let S^* → S_r. For both phase pressure drops scaled by the reference pressure drop, the result is linear trends with the inverse capillary number corrected for the phase viscosity times a saturation integral slope termed C_w and C_o for the respective phases. Expanding the capillary number shows that the relation is linear with inverse total velocity:$$\frac{\Delta p_w}{\Delta p_{w,\mathrm{ref}}}= 1+\frac{1}{N_0\left(\frac{{\mu }_w}{{\mu }_m}\right)}{C}_w=1+C_w\frac{{\sigma }_{\mathrm{ow}}}{{\mu }_{w}L\sqrt{\frac{\varphi }{K}}}\frac{1}{v_T},$$(76)

$$\frac{\Delta p_o}{\Delta p_{o,\mathrm{ref}}}= 1+\frac{1}{N_0\left(\frac{{\mu }_o}{{\mu }_m}\right)}{C}_o=1+C_o\frac{{\sigma }_{\mathrm{ow}}}{L\sqrt{\frac{\varphi }{K}}{\mu }_o}\frac{1}{v_T}.$$(77)

The saturation term slopes are defined as:$$C_w={\int }_{S=S_{\mathrm{eq}}}^{S_r}\left[k_{\mathrm{rw}}-k_{\mathrm{rw}}\left(S_{\mathrm{r}}\right)\right]\frac{\left(1-f_w\right)\frac{\mathrm{d}J}{\mathrm{d}S}}{F-f_w}\mathrm{d}S={\int }_{s_w=s_w^{\mathrm{eq}}}^{s_w^r}\left[k_{\mathrm{rw}}-k_{\mathrm{rw}}\left(s_w^r\right)\right]\frac{\left(1-f_w\right)\frac{\mathrm{d}J}{\mathrm{d}{s}_w}}{F-f_w}\mathrm{d}{s}_w,$$(78)

$$C_o={\int }_{S_{\mathrm{eq}}}^{S_r}\left[k_{\mathrm{ro}}-k_{\mathrm{ro}}\left(S_r\right)\right]\frac{f_w\frac{\mathrm{d}J}{\mathrm{d}S}}{F-f_w}\mathrm{d}S={\int }_{s_w=s_w^{\mathrm{eq}}}^{s_w^r}\left[k_{\mathrm{ro}}-k_{\mathrm{ro}}\left(S_r\right)\right]\frac{f_w\frac{\mathrm{d}J}{\mathrm{d}{s}_w}}{F-f_w}\mathrm{d}{s}_w.$$(79)

Equivalent to (68) we can use (76) and (77) to express the relations between effective and corrected relative permeabilities:$$\frac{1}{{\tilde{k}}_{\mathrm{rw}}}=\frac{1}{k_{\mathrm{rw}}\left(S_{\mathrm{r}}\right)}\left[ 1+\frac{1}{N_0\left(\frac{{\mu }_w}{{\mu }_m}\right)}{C}_w\right]=\frac{1}{k_{\mathrm{rw}}\left(S_{\mathrm{r}}\right)}\left[1+C_w\frac{{\sigma }_{\mathrm{ow}}}{L\sqrt{\frac{\varphi }{K}}{\mu }_w}\frac{1}{v_T}\right],$$(80)

$$\frac{1}{{\tilde{k}}_{\mathrm{ro}}}=\frac{1}{k_{\mathrm{ro}}\left(S_{\mathrm{r}}\right)}\left[ 1+\frac{1}{N_0\left(\frac{{\mu }_o}{{\mu }_m}\right)}{C}_o\right]=\frac{1}{k_{\mathrm{ro}}\left(S_{\mathrm{r}}\right)}\left[1+C_o\frac{{\sigma }_{\mathrm{ow}}}{L\sqrt{\frac{\varphi }{K}}{\mu }_o}\frac{1}{v_T}\right],$$(81)

which shows that if we plot the inverse of effective relative permeability against inverse capillary number or inverse velocity, we get a straight line trend with intercept at the ‘true’ relative permeability of each phase.

When calculating effective relative permeabilities for a given phase based on its own pressure drop data, we have shown that the intercept method follows where the inverse effective relative permeabilities plotted against inverse rate converge linearly to the correct inverse relative permeability without end effects. Assume now therefore that for a given phase the effective relative permeability is calculated using the pressure drop of the opposite phase. As shown in (71) this is directly related to that phase’s normalized pressure drop and combined with (76) and (77) we obtain:$$\frac{1}{{\tilde{k}}_{\mathrm{rw}}}=\frac{1}{k_{\mathrm{rw}}\left(S_{\mathrm{r}}\right)}\left[ \frac{\Delta p_o}{\Delta p_{\mathrm{ref}}}\right]=\frac{1}{k_{\mathrm{rw}}\left(S_{\mathrm{r}}\right)}\left[ 1+\frac{1}{N_0\left(\frac{{\mu }_o}{{\mu }_m}\right)}{C}_o\right],$$(82)

$$\frac{1}{{\tilde{k}}_{\mathrm{ro}}}=\frac{1}{k_{\mathrm{ro}}\left(S_{\mathrm{r}}\right)} \left[\frac{\Delta p_w}{\Delta p_{\mathrm{ref}}}\right]=\frac{1}{k_{\mathrm{ro}}\left(S_{\mathrm{r}}\right)}\left[ 1+\frac{1}{N_0\left(\frac{{\mu }_w}{{\mu }_m}\right)}{C}_w\right].$$(83)

We have thus demonstrated that the intercept method still holds: the inverse relative permeability is linear with the inverse capillary number and converges to the correct inverse relative permeability. However, the slope of the linear data will correspond to the phase with continuous pressure.

2.8 Single phase injection

An important special case is single phase injection (fractions F equal 0 or 1) which is used to determine critical saturations and relative permeability end points. As before, average saturation is determined by (50). The relative permeability of the phase not flowing cannot be determined directly, but will by definition be zero at steady state for the corrected saturation $s_w^r$. We now measure the pressure of the injected phase and its relative permeability is estimated using (80) if water is injected and (81) if oil is injected.

3 Results and discussion

Average (normalized) saturations can be calculated for any flow conditions based on (40) where the inlet saturation S₁ is determined from (42). Absolute average saturation then follows using (23). Pressure drops are calculated from (54), (55) and reference pressure drop from (64).

For mixed or unknown wettability, at a given flow fraction, the less-wetting phase can be identified as the phase whose saturation increases upon increased injection total rate. The measured pressure drop at that fraction is then associated with the less wetting phase. From our model we therefore choose to only report the non-wetting phase pressure drop.

Effective relative permeabilities are calculated based on the nonwetting phase pressure drop. These data are valid under any flow conditions and can be compared directly with steady state solutions from any commercial simulator. In our case the nonwetting phase is oil and we use (67) for effective oil relative permeability ${\tilde{k}}_{\mathrm{ro}}$ and (70) for effective water relative permeability ${\tilde{k}}_{\mathrm{rw}}$. An important exception is water flooding where the capillary pressure becomes zero or negative along the core and water will have higher pressure drop than the oil.

In addition, the analytical solutions corresponding to the intercept method are presented. They are valid only at high rates and will expectedly differ from the general solution at low rates. Selected effective relative permeability points are also calculated based on Sendra v2018.2.5 for validation.

3.1 Input data

Reference case input parameters are displayed in Table 1. As indicated, the oil-water flow parameters were based on measurements from Kleppe and Morse (1974) from a Berea sandstone core. Correlations with sufficient parameters to match the data were selected, such as the Andersen et al. (2017b) J-correlation for scaled capillary pressure and an extended Corey correlation (Brooks and Corey, 1964) for relative permeability with Corey exponents linearly dependent on saturation:$$J\left(S\right)=\frac{J_1}{{\left(1+k_{1}S\right)}^{n_1}}-\frac{J_2}{{\left(1+k_2\left(1-S\right)\right)}^{n_2}}+J_3,$$(84)

$$k_{\mathrm{rw}}=k_{\mathrm{rw}}^{\mathrm{*}}{\left(S\right)}^{n_w}, \hspace{1em}{k}_{\mathrm{rw}}=k_{\mathrm{ro}}^{\mathrm{*}}{\left(1-S\right)}^{n_o},$$(85)

$$n_w=n_{w1}S+n_{w2}\left(1-S\right),\hspace{1em} n_o=n_{o1}S+n_{o2}\left(1-S\right).$$(86)

All the functions were expressed using the normalized saturation S from (23). The J-function derivative, which is central in the solution is easily evaluated explicitly as:$$\frac{\mathrm{d}J}{\mathrm{d}S} =-\frac{J_1{n}_1{k}_1}{{\left(1+k_{1}S\right)}^{n_1+1}}-\frac{J_2{n}_2{k}_2}{{\left(1+k_2\left(1-S\right)\right)}^{n_2+1}}.$$(87)

The relative permeabilities and scaled and unscaled capillary pressure can be seen in Figure 2 comparing correlations with the experimental data. Fractional flow functions with three choices of oil viscosity are shown in Figure 3.

3.2 Validation

The analytical solution was validated by comparison of results with full transient numerical simulations from the commercial software Sendra v2018.2.5. This program solves the flow equations fully implicit. The core was discretized using 500 grid cells to capture steep saturation gradients. At each injection condition the program was run until a visibly flat production profile was observed, indicating steady state.

Using the base case parameters, two flow fractions are considered: F = 0.50 and F = 0.05 and injection rate is varied for both cases using 0.1, 1, 10 and 100 pore volumes per day (PV/d). The resulting steady state saturation profiles are shown in Figure 4 comparing the analytical and numerical solutions. For comparison, the saturation integrals in the derived methods were calculated using 10 000 saturations in the interval between S_r and S_eq, sufficient for the integrals to converge.

For all rates and flow fractions the saturation distribution solutions from the analytical solution (28) and the commercial software (Sendra) overlap, verifying our procedure. All the saturation distributions approach the saturation $s_w^{\mathrm{eq}}=0.61$ at Y = 0. For a given flow fraction the saturation distributions approach the saturation $s_w^r$ such that f_w$\left(s_w^r\right)$ = F. Considering Figure 3 we see that when F = 0.5 and F = 0.05, then $s_w^r\approx 0.545$ and $s_w^r$ ≈ 0.475, respectively, which is consistent with the profiles in Figure 4. For high rates the distributions stabilize at $s_w^r$ within the core a certain distance from the outlet. For low rates (such as 0.1 PV/d in these examples), however, the capillary forces dominate and shift the saturations towards $s_w^{\mathrm{eq}}$ and no saturations in the core are close to $s_w^r$. Direct application of Darcy’s law to calculate relative permeability at steady state assumes $sw\left(Y\right)\mathrm{ }=\mathrm{ }{s}_w^r$ which we see is only reasonable at very high rates.

3.3 Profile shapes

The analytical solution predicts that variation of N₀ only compresses or expands the saturation profile as long as the saturation functions, viscosity ratio and flow fraction are not changed, see (44). This is demonstrated by varying the core length L, mean viscosity μ_m, total velocity v_T and permeability K (which all are directly involved in N₀, see (27)) around the base case with F = 0.05 and rate 1 PV/d to give same change in N₀ (an increase of N₀ by a factor 5). The corresponding saturation profiles are shown in Figure 5a together with the reference case. As expected, all the parameter variations increasing N₀ by the factor 5 result in the same scaled profile which is only a compression of the base profile. For example, the saturation s_w = 0.5 has position Y ≈ 0.9 in the base case and $Y\approx \frac{0.9}{5}$ = 0.18 in the new cases.

In comparison, in Figure 5b we see that variations in individual viscosities (thus changing the viscosity ratio) and flow fraction affect both magnitude and shape of the profiles as they do not overlap upon scaling. In particular, different values of the end point $s_w^r$ are obtained in each case.

3.4 Determining the n-parameter for approximate profiles

The n parameter appears in the scaling number Y_cee. To find a good choice for n we want to see that the actual numerical profiles are well approximated by the polynomial function. This was considered by varying the injected fraction F (with values 0.01, 0.05 and 0.5) and the oil viscosity μ_o (with values 1, 5 and 25 times the reference value of 2.3 cP) to get distinct profile shapes. The normalized saturation profiles S(Y) (blue curves) were plotted together with the corresponding approximated profiles (orange curves) in Figure 6 for n = 5 which was found to be a good choice. The profiles are plotted over the range Y = 0 to 10Y_av. The positions Y_av and Y_cee = (n + 1)Y_av are marked with vertical dashed lines (left and right, respectively). Also a dashed/dotted line indicating the normalized saturation S_r (obtained at infinite distance) is shown. The main observations are:

• The major part of the saturation profile is located behind Y_cee including 90–95% of the saturations closest to S_eq.

• Saturations located at Y > Y_cee are very close to S_r.

• We can state end effects to be severe when Y_cee > 1 since then all saturations deviate from S_r more than 5% the magnitude of the end effect saturation interval |S_r − S_eq|.

In terms of the impact the fraction and viscosities have, we see that when the fraction F is lowered, S_r decreases farther from S_eq. This expands the saturation range and the extent of end effects. We also see that increasing the oil viscosity alters the mobility ratio to reduce S_r and expands the saturation range. Although this mechanism increases the end effect, the increased viscous forces work to compress the profile. At low F the increase in oil viscosity has more impact on compressing the end effect profile.

3.5 Intercept method

3.5.1 Demonstration of analytical solution/intercept method

Consider the reference case for the four combinations of flow fractions F = 0.05 and 0.01 and oil viscosities 2.3 cP and 5 × 2.3 cP. Average saturation and inverse effective oil and water relative permeabilities, $\frac{1}{{\tilde{k}}_{\mathrm{rw}}}$ and $\frac{1}{{\tilde{k}}_{\mathrm{ro}}}$, are plotted against inverse velocity $\frac{1}{v_T}$ (Fig. 7). Flow rates are used in the range 5 PV/d to 1000 PV/d. The following is seen:

• For all cases, once the rates are sufficiently high, the results overlap completely with the linear analytical solution. As $\frac{1}{v_T}\to 0$, the solutions reach the correct saturation $s_w^r$ and relative permeability points $\frac{1}{k_{\mathrm{ri}}\left(s_w^r\right)}$ without end effects, as marked on the y-axis with large symbols.

• At low rates the two solutions do not overlap. The analytical solution is only valid when the end effect profile has not reached the outlet significantly which is less likely at low rates.

• When plotting the relative permeability data on two scales such that two points overlap, all the points overlap, i.e. the profile shapes are identical.

• Since the trends are linear, and the linear trends appear valid over a wide range of rates, only a few rates are necessary to determine the corrected values.

• The above statements were validated for average saturation and relative permeability by simulations with the software Sendra for the case F = 0.01, μ_o = 2.3 cP.

3.5.2 Scaling of end effects with dimensionless numbers

We have seen that expressing average saturation and inverse relative permeabilities as function of inverse velocity results in linear relations. The velocity dependent terms directly express the impact of end effects on those properties and are reasonable to apply as dimensionless numbers for scaling the behavior.

In Figure 8a we plot ${\overline{s}}_w-s_w^r$ against the number $\frac{C_s}{N_0}$ based on different flow fractions (F = 0.9–0.005) and injection rates and in Figure 8c we do the same varying oil viscosities (from 5⁰ to 5⁴ times the reference value) and injection rates. For a given flow fraction, C_s is constant while N₀ varies with rate. Oil viscosity changes both terms.

At high rates ${\overline{s}}_w-s_w^r$ is identical to $\frac{C_s}{N_0}$, which follows directly from the analytical solution, see (50), and all cases fall on this straight line. The dimensionless number $\frac{C_s}{N_0}$ thus scales saturation end effects and accounts for all relevant input parameters. For low rates the cases deviate from the straight line with different trends. The number $\frac{C_s}{N_0}$ directly expresses the change in average saturation due to end effects. As end effects also can reduce the average saturation, the straight line can continue to negative values. Values of $\frac{C_s}{N_0}$ can be considered ‘large’ if it changes the average saturation by 0.01 or more (for comparison the mobile saturation range in our data is 0.31) and negligible for $\frac{C_s}{N_0}<0.001$.

We also plot $\frac{k_{\mathrm{ri}}\left(s_w^r\right)}{{\tilde{k}}_{\mathrm{ri}}}-1$ against $\frac{C_o}{N_0\left(\frac{{\mu }_o}{{\mu }_m}\right)}$ in Figure 8b for different flow fractions and Figure 8d for different oil viscosities. From (80) to (83) we see that the two expressions are identical and form a straight line from zero at high rates (low $\frac{C_o}{N_0\left(\frac{{\mu }_o}{{\mu }_m}\right)}$). At low rates the data deviate from the straight line.

On the straight line, the term $\frac{C_o}{N_0\left(\frac{{\mu }_o}{{\mu }_m}\right)}$ expresses how much the effective relative permeability is reduced compared to the corrected relative permeability as given by ${\tilde{k}}_{\mathrm{ri}}=\frac{1}{1+\left[\frac{C_o}{N_0\left(\frac{{\mu }_o}{{\mu }_m}\right)}\right]}{k}_{\mathrm{ri}}\left(s_w^r\right)$. If $\frac{C_o}{N_0\left(\frac{{\mu }_o}{{\mu }_m}\right)}=0$, they are identical, while if nonzero they differ. A value of $\frac{C_o}{N_0\left(\frac{{\mu }_o}{{\mu }_m}\right)}$ greater than 1 can be considered large since then ${\tilde{k}}_{\mathrm{ri}}<\frac{1}{2}{k}_{\mathrm{ri}}\left(s_w^r\right)$. Since the measured pressure drop is in the less wetting phase which has the highest pressure drop, $\frac{C_o}{N_0\left(\frac{{\mu }_o}{{\mu }_m}\right)}$ will be positive and the effective relative permeabilities will be lower than the corrected relative permeabilities, i.e. ${\tilde{k}}_{\mathrm{ri}}<k_{\mathrm{ri}}\left(s_w^r\right)$. This is seen in the graph as $\frac{k_{\mathrm{ri}}\left(s_w^r\right)}{{\tilde{k}}_{\mathrm{ri}}}-1>0$.

If conditions were such that water was non-wetting at a given fraction (in a mixed-wet or oil-wet system) then we should modify the plot to show $\frac{k_{\mathrm{ri}}\left(s_w^r\right)}{{\tilde{k}}_{\mathrm{ri}}}-1$ against $\frac{C_w}{N_0\left(\frac{{\mu }_w}{{\mu }_m}\right)}$.

An interesting observation is that when the curves deviate from the straight line, they all fall below it, indicating that the impact is less than if the entire end effect profile was within the core. That is reasonable since the profile that is not within the core does not contribute to change the observations, and since at lower rates more of the profile is outside the core, the added impact on the saturation and relative permeability becomes less.

Note also that the dimensionless number that scales saturation, $\frac{C_s}{N_0}$, is different from the one scaling relative permeability, $\frac{C_w}{N_0\left(\frac{{\mu }_w}{{\mu }_m}\right)}$ or $\frac{C_o}{N_0\left(\frac{{\mu }_o}{{\mu }_m}\right)}$. As mentioned, the dimensionless number which is used for scaling relative permeability is also determined by which phase is more wetting at the given fraction (i.e. whether $s_w^r$ is above or below $s_w^{\mathrm{eq}}$ which further depends on wettability, flow fraction and viscosity ratio).

To see more systematically when the solutions deviate from the analytical solution (the linear model) we plot ${\overline{s}}_w^{\mathrm{AnSol}}-{\overline{s}}_w$ and $\frac{{\left(\frac{k_{\mathrm{ri}}\left(s_w^r\right)}{{\tilde{k}}_{\mathrm{ri}}}\right)}^{\mathrm{AnSol}}}{\left(\frac{k_{\mathrm{ri}}\left(s_w^r\right)}{{\tilde{k}}_{\mathrm{ri}}}\right)}=\frac{{\tilde{k}}_{\mathrm{ri}}}{{\tilde{k}}_{\mathrm{ri}}^{\mathrm{AnSol}}}$ for the aforementioned cases in Figure 9, using the same dimensionless numbers $\frac{C_s}{N_0}$ and $\frac{C_o}{N_0\left(\frac{{\mu }_o}{{\mu }_m}\right)}$, respectively. When the analytical solution is valid, these terms are 0 and 1, respectively. This is the case at low values of the respective dimensionless numbers. The onset of deviation in saturation becomes visible for 0.005 < $\frac{C_s}{N_0}$ < 0.02 for the cases with variation in F and 0.008 < $\frac{C_s}{N_0}$ < 0.03 for the cases with variation in oil viscosity. The onset of deviation in relative permeability trends $\frac{{\tilde{k}}_{\mathrm{ri}}}{{\tilde{k}}_{\mathrm{ri}}^{\mathrm{AnSol}}}$ from 1 is visible for 0.8 < $\frac{C_o}{N_0\left(\frac{{\mu }_o}{{\mu }_m}\right)}$ < 5 in the cases with variation in F and 0.8 < $\frac{C_o}{N_0\left(\frac{{\mu }_o}{{\mu }_m}\right)}$ < 3 for the cases with variation in oil viscosity. The onset thus occurs at values of the dimensionless numbers varying by an order of magnitude. Further, both dimensionless numbers appear directly linearly in their separate relations (for saturation or relative permeability). However, their magnitudes at which we observe deviation from the linear relations are very different, spanning three orders of magnitude. Hence, although the dimensionless numbers are excellent at scaling the end effects, they do not accurately determine when the linear analytical solutions cease to be valid.

However, we found that Y_cee is a good measure of when significant changes in saturation have occurred along the entire core and previously found an appropriate value of n = 5 to apply in this number from matching saturation profiles, see Figure 6. We now present the above results of ${\overline{s}}_w^{\mathrm{AnSol}}-{\overline{s}}_w$ and $\frac{{\tilde{k}}_{\mathrm{ri}}}{{\tilde{k}}_{\mathrm{ri}}^{\mathrm{AnSol}}}$ plotted against Y_cee in Figure 10. For all cases and both the saturation and relative permeability data we see that the onset of deviation starts at Y_cee ≈ 1 with very little variation (a factor 1.2, i.e. 0.08 orders of magnitude). As we have reduced the span from 1 order of magnitude for individual terms and from 3 orders of magnitude considering both terms ${\overline{s}}_w^{\mathrm{AnSol}}-{\overline{s}}_w$ and $\frac{{\tilde{k}}_{\mathrm{ri}}}{{\tilde{k}}_{\mathrm{ri}}^{\mathrm{AnSol}}}$ to 0.08 orders we can claim that Y_cee < 1 is a reliable criterion for when the intercept method is valid.

3.6 Appearance of rate-dependent relative permeabilities

End effects can affect both the steady state average saturation and effective relative permeability points, which deviate from the true relative permeability functions. Effectively this can result in an appearance of rate-dependent relative permeability functions meaning that if they are measured at a fixed rate, the rate will affect which curves are obtained and not offer information about their correctness.

Based on the reference case relative permeability functions and capillary pressure as input, effective relative permeabilities were calculated and plotted against average saturation. For a given injection total rate, 40 evenly spaced fractions were used between F = 0.005 and 0.995, each producing a relative permeability point. This was performed for nine injection total rates ranging from 0.1 to 1000 PV/d varied by the same factor 10^0.5. Since not only rate, but the dimensionless number N₀ determines the impact of end effects, the effective relative permeability curves were labeled by the corresponding value of N₀ and are shown in Figure 11. High rate corresponds to high N₀ and opposite. The input relative permeabilities, which are the true curves of the system, are included for comparison with the effective curves.

Notably, at low N₀ (low total rates) the average saturations are shifted towards $s_w^{\mathrm{eq}}$ = 0.61, leading to narrower saturation ranges. Each calculated relative permeability point is reduced by the high pressure drop in the nonwetting phase. At high N₀ (high total rates) the relative permeability curves converge and become independent of injection total rate (and N₀). They then also agree with the input curves (“true rel perms”) used to generate the data.

We do not argue that relative permeabilities physically cannot exhibit rate-dependence. For example, micromodel observations indicate the possible existence of a flow regime dependent of the conventional capillary number (velocity times viscosity divided by interfacial tension) (Valavanides, 2018).

3.7 Calculation of relative permeability and capillary pressure from experiments

We consider steady state experiment data where two phases are injected at different flow fractions F^j and at each fraction different total velocities $v_T^k$ are applied. Potentially such tests can be applied using different cores and even fluids, under the assumption that the saturation functions remain the same. For example, different core tests can be implemented at constant total velocity using the same schedule of injected fractions (Henderson et al., 1998; Virnovsky et al., 1998). The cores used can have different length L, porosity ϕ, permeability K. Use of different fluids (or just measuring a different value for the same fluids in a different test) can change viscosities μ_i and interfacial tension σ_ow. Average saturation ${\overline{s}}_w^{j,k}$ and pressure drop Δp^j,k are measured at each fraction and velocity.

3.7.1 Simple procedure: Calculation of only relative permeability with constant fluid and rock properties

1. Convert the pressure drop data to effective relative permeability measurements ${\tilde{k}}_{\mathrm{ri}}^{j,k}$.

$${\tilde{k}}_{\mathrm{ri}}^{j,k}=\frac{\mathrm{L\varphi }{\mu }_i}{K}\frac{F_i^j{v}_T^{j,k}}{\Delta p^{j,k}}.$$(88)

1. For each fraction F^j plot average saturation ${\overline{s}}_w^j$ against inverse total rate $\frac{1}{v_T}$ and inverse effective relative permeabilities $\frac{1}{{\tilde{k}}_{\mathrm{ri}}^{j,k}}$ against inverse total rate $\frac{1}{v_T}$. Determine the lines through the data:

$${\overline{s}}_w^j=s_w^{r,j}+m_s^j\frac{1}{v_T},$$(89)

$$\frac{1}{{\tilde{k}}_{\mathrm{ri}}^j}=\frac{1}{k_{\mathrm{ri}}\left(s_w^{r,j}\right)}+m_i^j\frac{1}{v_T}.$$(90)

Report the intercepts $s_w^{r,j},\frac{1}{k_{\mathrm{ri}}\left(s_w^{r,j}\right)}$. The slopes $m_s^j,m_i^j$ are not used.

1. The pairs $\left(s_w^{r,j},k_{\mathrm{ri}}\left(s_w^{r,j}\right)\right)$ for i = o, w are accurate relative permeability points. Match suitable relative permeability correlations (and critical saturations) to fit them.

3.7.2 Advanced procedure: Calculation of both relative permeability and capillary pressure from tests with distinct properties

1. For each fraction F^j and total velocity $v_T^k$, calculate the number $N_0^{j,k}={\left[\frac{v_{T}L\sqrt{\frac{\varphi }{K}}{\mu }_m}{{\sigma }_{\mathrm{ow}}}\right]}^{j,k}$ and the viscosity ratios ${\left(\frac{{\mu }_i}{{\mu }_m}\right)}^{j,k}$. The involved parameters are standard.

2. For each fraction F^j plot average saturation points ${\overline{s}}_w^{j,k}$ against $\frac{1}{N_0}$. Determine the line through the data (all or those at lowest $\frac{1}{N_0}$) by regression and report the slope $C_s^j$ and intercept $s_w^{r,j}$ as:

$${\overline{s}}_w^j=s_w^{r,j}+C_s^j\frac{1}{N_0}$$(91)

1. For each fraction F^j determine which phase is more wetting. If average water saturation decreases at lower $\frac{1}{N_0}$, water is wetting at that fraction. If opposite, oil is wetting. The measured pressure drop corresponds to the less wetting phase.

2. Convert the pressure drop data to effective relative permeability measurements ${\tilde{k}}_{\mathrm{ri}}^{j,k}$:

$${\tilde{k}}_{\mathrm{ri}}^{j,k}={\left[\frac{\mathrm{L\varphi }{v}_T{\mu }_i{F}_i}{K\Delta p}\right]}^{j,k}.$$(92)

1. Assume first the fluids have constant viscosities. For each fraction F^j and rate $v_T^k$ plot inverse effective relative permeabilities $\frac{1}{{\tilde{k}}_{\mathrm{ri}}^{j,k}}$ against $\frac{1}{N_0}$. Determine the line through the points (all or those at lowest $\frac{1}{N_0}$) and report the slope $m_i^j$ and intercept $\frac{1}{k_{\mathrm{ri}}\left(s_w^{r,j}\right)}$ as:

$$\frac{1}{{\tilde{k}}_{\mathrm{ri}}^j}=\frac{1}{k_{\mathrm{ri}}\left(s_w^{r,j}\right)}\left[1+m_i^j\frac{1}{N_0}\right].$$(93)

Both relative permeability lines will have the same slope $m_o^j=m_w^j=m^j$. Depending on which phase is wetting at the given fraction we determine either $C_o^j$ or $C_w^j$:$$C_o^j=m^j{\left(\frac{{\mu }_o}{{\mu }_m}\right)}^j,\hspace{1em}\left(\mathrm{if} \mathrm{oil} \mathrm{less} \mathrm{wetting}\right),$$(94)

$$C_w^j=m^j{\left(\frac{{\mu }_w}{{\mu }_m}\right)}^j,\hspace{1em}\left(\mathrm{if} \mathrm{water} \mathrm{less} \mathrm{wetting}\right).$$(95)

If the viscosities vary between tests, the data must be plotted against $\frac{1}{N_0\left(\frac{{\mu }_w}{{\mu }_m}\right)}$ for the fractions where water is less wetting and against $\frac{1}{N_0\left(\frac{{\mu }_o}{{\mu }_m}\right)}$ where oil is less wetting.

1. The pairs $\left(s_w^{r,j},k_{\mathrm{ri}}\left(s_w^{r,j}\right)\right)$ for i = o, w are accurately determined relative permeability points.

2. Match suitable relative permeability correlations (and critical saturations) to the relative permeability points.

3. The slopes $C_s^j,C_o^j,C_w^j$ indicate how strongly capillary forces affect saturation and pressure drop at the relevant saturation ranges $s_w^{\mathrm{eq}}$ to $s_w^{r,j}$ and constrain the choice of capillary pressure functions. Their values only depend on the saturation functions and are independent of fluid and core properties.

4. Assume a capillary pressure correlation with tuning parameters. Based on the calculated relative permeability correlations, select tuning parameters and calculate slope constants $C_s^j,C_o^j,C_w^j$. The tuning parameters giving least relative error between predicted and measured values result in the capillary pressure J-function best explaining the data.

3.7.3 Interpretation of experimental data I

Virnovsky et al. (1998) conducted three oil-water steady state tests on a Berea sandstone core. The core was initially fully water saturated and oil fractions were increased. A fixed total rate was used in each test of 0.2, 0.5 or 5 cc/min. The absolute permeability changed between the three tests, hence we used N₀ to account for the differences. From the three velocities at each fraction we follow the ‘advanced’ procedure. Average saturations and inverse effective relative permeabilities were plotted against $\frac{1}{N_0}$ in Figures 12a–12d, respectively to find the intercepts $\left(s_w^{r,j},k_{\mathrm{ri}}\left(s_w^{r,j}\right)\right)$ and the slope constants C_s, C_o, C_w which are listed in Table 2. Linear trends in the points were clearly visible for most cases.

For the 5 highest fractions F water saturation increased with rate, indicating capillary forces trapped oil then, while for the lowest 3 fractions water saturation decreased with increased rate, indicated capillary forces then trapped water. Equivalently C_s, the slope of average saturation vs. $\frac{1}{N_0}$, was negative for the highest fractions and positive for the lowest fractions. Although there was some uncertainty in the sign for some central fractions it seems clear that $s_w^{\mathrm{eq}}$ was located in the intermediate saturation range.

As described in the theory, we can scale end effects by plotting ${\overline{s}}_w-s_w^r$ against $\frac{C_s}{N_0}$ and $\frac{k_{\mathrm{ri}}\left(s_w^r\right)}{{\tilde{k}}_{\mathrm{ri}}}-1$ against $\frac{C_o}{N_0\left(\frac{{\mu }_o}{{\mu }_m}\right)}$ (when oil is less wetting) and against $\frac{C_w}{N_0\left(\frac{{\mu }_w}{{\mu }_m}\right)}$ (when water is less wetting). This was done for all the data points and is shown in Figure 13. Average saturation, water relative permeability and oil relative permeability fall on a straight line in their respective plot.

The relative permeability points were fitted with the correlation (85) and (86), see Figure 14a. Next, the J-function was determined, see Figure 14b, using the correlation (84) by optimizing the match of slope parameters: The C_s values were well matched, see Figure 15a, while it was more difficult to match the C_o, C_w parameters which at optimum were roughly 2–10 times lower than those derived from the measurements, see Figure 15b. This means that the model with the tuned curves would capture well what the average saturation and effective relative permeability values would approach at each fraction, and that the slope of average saturation with inverse total rate would vary as observed experimentally. However, the effective relative permeabilities would change less in the model than observed experimentally.

Near identical J-functions resulted from different optimizations, indicating that the data and approach give consistent results. The optimal curve parameters and the system parameters are found in Table 3. Comparison with a measured drainage capillary pressure curve indicated similar magnitude and shape, however negative capillary pressures at high saturations resulted from the model, required to explain the trends in saturation with rate. Also, a higher residual water saturation was found which seems more consistent with the flooding data.

The model can under any conditions only be sensitive to saturations on the interval between $s_w^{\mathrm{eq}}$ and $s_w^r$. As $s_w^{\mathrm{eq}}$ was located centrally, the extreme $s_w^r$ values corresponding to the lowest and highest fractions gave the total saturation range where we have accurate information. As seen in Table 2 and Figure 14, this is between saturations 0.326 and 0.967 which is practically the entire mobile saturation range.

3.7.4 Interpretation of experimental data II

Henderson et al. (1998) measured relative permeability to gas and condensate at 7 flow fractions between 0.048 and 0.29 (corresponding to condensate-gas-ratio 0.05–0.4) and 4 total velocities on Berea sandstone with presence of connate water saturation. Condensate c and gas g were treated as ‘oil’ and ‘water’ in the above theory. Plotting average condensate saturations ${\overline{s}}_c$ and inverse effective relative permeabilities $\frac{1}{{\tilde{k}}_{\mathrm{rc}}}$ and $\frac{1}{{\tilde{k}}_{\mathrm{rg}}}$ against $\frac{1}{N_0}$ in Figure 16 indicates that the points at the three highest rates (lowest $\frac{1}{N_0}$) aligned well on straight lines. This was used to obtain the intercepts and resulting points $\left(s_c^r,k_{\mathrm{ri}}\left(s_c^r\right)\right)$ for i = c, g, as well as the slope parameters C_s, C_c, C_g. All these values are listed in Table 4 for their respective flow fractions.

For all the data, the points at the lowest rate (highest $\frac{1}{N_0}$) that did not fall on the line, had less impact on the saturation or relative permeability than if they had been on the line (the trends flattened). This is consistent with what was shown in Figure 8 where the end effects’ added impact lessens as more of the profile is outside the core.

Also, at each fraction, increased velocity at a fixed fraction increased the condensate saturation, or equivalently condensate saturation decreased with higher $\frac{1}{N_0}$, quantified by negative C_s. This indicated that capillary forces trapped gas under all conditions and gas wetted the rock more than condensate for the considered fractions.

Scaling the saturation and relative permeability points is shown in Figure 17 where we again see that all the high rate points fall well on the same straight line when plotted against the relevant dimensionless numbers. The low rate results fall between the line and the horizontal axis indicating that the impact is less than if they had continued on the line.

The intercept saturation and relative permeability points were collected and matched with correlations, as shown in Figure 18a. The slope parameters were then matched by tuning the J-function, seen in Figure 18b. Only the saturation range where the curves are reliable is shown. That includes the range between the lowest and highest $s_c^r$ (0.397–0.507), but also down to $s_c^{\mathrm{eq}}$ = 0.344 since all saturation profiles end at that saturation. The optimal match between the slope parameters from the tuned model and the values from the observations is shown in Figure 19. Both saturation and relative permeability slope parameters are well matched, with almost all values from the tuned model being less than 30% different from the values from the measurements. System parameters and correlation parameters from matching the data are listed in Table 5.

Henderson et al. (1998) interpreted their data as capillary number dependent relative permeabilities rather than being affected by capillary end effects. They suggested that the pressure drop across the core was orders of magnitude higher than the capillary pressure that could cause end effects. Although we cannot conclude that end effects were present or not, they explain the observations very consistently in that the data fall on the predicted straight lines at high rates, fall on the expected side of the line at low rates and we can obtain saturation functions that consistently explain all the data. Gupta and Maloney (2016) interpreted the same dataset to find corrected relative permeabilities.

4 Conclusion

A theory is derived from fundamental assumptions describing how capillary end effects affect average saturation and relative permeability calculation in steady state experiments. The method is valid for all saturations functions and wetting states.

We derive an intuitive and general ‘intercept’ method for correcting steady state relative permeability measurements for capillary end effects: plotting average saturation and inverse relative permeability against inverse rate will give linear trends at high rates and result in corrected relative permeability points when extrapolated to zero inverse rate (infinite rate). This is a formal proof and generalization of the method proposed by Gupta and Maloney (2016).

It is shown how the slopes of the lines are related to the saturation functions allowing to scale all data for all conditions to the same straight lines. Two dimensionless numbers are obtained that directly express when end effects become important in terms of (a) how much the average saturation is changed: C_s/N₀ and (b) how much effective relative permeabilities are reduced $\frac{C_i}{N_0}$.

A third dimensionless number Y_cee (expressing the scaled front position of the end effect region) is derived directly stating that the end effect profile has reached the inlet of the core significantly when Y_cee = 1. This number acts as a universal criterion for when the linear (intercept method) behavior begins. The intercept method is valid when Y_cee < 1.

All the dimensionless numbers contain a part depending only on saturation functions, injected flow fraction and viscosity ratio and a second part N₀ containing constant known fluid and rock parameters such as core length, porosity, interfacial tension, etc. The former parameters determine the saturation range and shape of the saturation profile, while the number N₀ determines how much the profile is compressed towards the outlet.

End effects cause the average saturations to be shifted towards $s_w^{\mathrm{eq}}$ (the saturation at which capillary pressure is zero) and the calculated relative permeabilities from pressure drop and Darcy’s law to be reduced in magnitude compared to the true relative permeabilities. Since the shift of the effective curves depends on the viscous vs. capillary forces, using different rates will result in different effective curves when capillary end effects are significant. This gives a false impression of rate-dependent relative permeabilities.

Methodologies for deriving relative permeability and capillary pressure systematically and consistently, even based on combining data from tests with different fluid and core properties, are presented and demonstrated on two datasets from the literature.

Even with access to only one phase pressure drop of two injected phases it is shown that the intercept method holds for both phase relative permeabilities. This is the standard regarding what pressure data is measured.

Acknowledgments

The author acknowledges the Research Council of Norway and the industry partners, ConocoPhillips Skandinavia AS, Aker BP ASA, Vår Energi AS, Equinor ASA, Neptune Energy Norge AS, Lundin Norway AS, Halliburton AS, Schlumberger Norge AS, and Wintershall DEA, of The National IOR Centre of Norway for support.

References

All Tables

All Figures