Steady state relative permeability experiments with capillary end effects: Analytical solutions including derivation of the intercept method

Steady state relative permeability experiments are performed by coinjection of two (or more) immiscible fluids. The relative permeabilities can be calculated directly from the stabilized pressure drop and saturation of the core if capillary end effects and transient effects are negligible. In most cases such conditions are difficult to obtain. This work presents an analytical solution in form of explicit expressions for the spatial profiles of pressure gradients and saturation, average saturation and pressure drop for a core being injected simultaneously with two phases at steady state when capillary end effects are significant. When arbitrary saturation functions are applied, such parameters and distributions can only by obtained by numerical integration. By assumption of simplified saturation function correlations the differential equation describing steady state can be integrated. A new dimensionless capillary number N is obtained which contains the fluid and rock parameters, but also the saturation function parameters (relative permeability and capillary pressure), fluid viscosities, injected flow fraction, total flow rate and more. It is shown that when this number is of magnitude 1, end effects cover parts of the core, but parts of the core are also unaffected. For N > 10 the end effects are negligible, while for N < 10 end effects are dominant. This paper gives the first formal proof of the intercept method from basic assumptions. It is shown that when the inlet saturation is sufficiently close to that of a no capillary pressure situation; the average saturation changes linearly with the inverse of total rate towards the saturation corresponding to no capillary forces; also, the pressure drop divided by the pressure drop of a no end effect situation goes linearly towards 1 with the inverse of total rate.


Introduction
Relative permeabilities are parameters describing how the mobility of a phase is affected in presence of other fluids in a porous medium. Accurate measurements of relative permeability are required for making reliable predictions and decisions at field scale. Such measurements are traditionally performed either by the unsteady state method where one fluid is injected to displace the other or the steady state method where both fluids are co-injected. The former method is similar to the displacement taking place in the reservoir, while stable and uniform flow can be achieved in the latter. The steady state method will be the focus of this paper.
Core flooding where several phases are involved will be affected by capillary pressure. Leverett (1941) stated that the natural outlet boundary condition is a zero capillary pressure. This follows from the phase pressures being continuous towards the outlet where both phases are produced and that the radius of curvature goes from a small value in the pore space to (practically) infinity. This condition forces the saturation at the outlet to be fixed as defined by the capillary pressure curve (the imbibition curve if it is an imbibition process) and can be observed as a fluid accumulation near the outlet (Richardson et al., 1952). The amount of accumulation will depend on wettability, injected water fraction, advective forces and capillary forces. Increasing the injection rate will in such cases affect the average saturation for the same injected fluid fraction. Especially, for water injection this affects estimation of the critical oil saturation and end point water relative permeability. Apparently rate-dependent relative permeabilities do however tend to give consistent results once the end effects have been corrected for (Osoba et al., 1951;Chen and Wood, 2001). Virnovsky et al. (1995) developed a model which expressed the relative permeabilities and capillary pressure as functions of the responses of average saturation and pressure drop to changes in total rate. To make full use of their model they would require measurements for many different rates at a given fraction. Huang and Honarpour (1998) developed a model for oilflooding to displace water in a strongly water-wet system and obtained analytical solutions by means of Corey-Burdine equations relating the capillary pressure and relative permeability functions. Although analytical expressions could be defined for saturation profiles and average saturation, they depended on the inlet saturation which needed calculation from an algebraic equation. Gupta and Maloney (2016) recently developed the intercept method. This method argues that stabilized pressure drop and average saturation measurements during steady state tests vary systematically with rate; if end effects are limited to within the core, then the average saturation s w and pressure drop divided by rate Δp=Q plotted against inverse rate 1 Q will linearly approach constant values s rep w and � Δp Q � rep at inverse rate equal zero (with slopes a and b), representative of a system without end effects.
This is especially important since measurements taken within the limits of experimental metering, sample integrity and Darcy's law could be extrapolated to the correct result. However, in their derivation they mainly assumed that the role of the capillary pressure was a fixed pressure drop and did not derive their method for arbitrary saturation functions. For a review of the intercept method and its applications, see Reed and Maas (2018). Andersen et al. (2017) assumed saturation functions of a Corey type for water relative permeability and scaled capillary pressure and derived explicit analytical solutions for oil displacement by waterflooding affected by end effects in a mixed-wet system. A dimensionless capillary number incorporating the saturation function parameters was derived such that for its critical value of 1 the end effect region exactly reached across the core from outlet to inlet. For higher values of this number (i.e. high rates) the intercept method was obtained analytically for this system (waterflooding). In addition, the method explained how saturation and pressure drop varied when the theoretical end effect region also went beyond the core length. The model was used to interpret experimental data in Andersen et al. (2020) and estimate the full saturation functions and not only one point of the relative permeabilities as originally intended.
In this work we aim to derive analytical solutions for steady state relative permeability measurements in form of explicit expressions. This requires the use of specific and limited correlations forms, but allows us to derive illustrative and theoretical results. By taking a simplest possible approach with linear saturation functions we underline that the correlations are not considered flexible enough to represent experimental data, but capture physical behavior related to saturation change on mobility and capillary pressure and how this interplays with the core and fluid properties and typical experimental control settings such as total injection rate and fluid fractions. We also show that the intercept method is obtained from this solution.
The paper is structured as follows: a) We present equations for steady state flow during co-injection of two phases. Assuming appropriate forms of the saturation functions we derive explicit analytical solutions for spatial saturation and pressure gradient profiles and average saturation and phase pressure drops. b) Illustrative and theoretical results of the analytical solution are presented. Particularly the solution demonstrates a derivation of the intercept method. c) The paper is summarized by conclusions.

General model description
The mathematical description of 1D incompressible and immiscible flow of oil (o) and water (w) in a porous homogeneous medium is given by: where φ is porosity, s i saturation, 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: The total Darcy velocity u T is defined as: It follows from adding the transport equations in (2) that: The water phase equation (2) can then be expressed with variables u T ; s w as: where f w ¼ λw λwþ λo is the fractional flow function.

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) with a total Darcy flux u T : The flux of a given phase is composed of both an advective and capillary component. Hence, we note that F does not correspond to f w ðx ¼ 0Þ unless capillary forces can be ignored. From (8) we write this boundary condition as: The outlet boundary condition is described by a zero capillary pressure: The initial condition is considered the steady state of a previous injected fraction, F pre . If F pre ¼ 0 this state corresponds to initial water saturation. We assume an imbibition process is considered where water saturation will increase with time as implemented by setting F > F pre , i.e. injecting a higher fraction of water.

Steady state
At steady state we have no changes with time in the system, i.e.: 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, saturations and pressures will be taken as function of spatial coordinate alone, e.g. s w ¼ s w ðxÞ. (8) can be written as: 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 gives: This is equivalent to: Using that d x P c ¼ P ' c ðs w Þd x s w , we can solve (14) with respect to the saturation gradient: 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 and dimensionless Leverett J-function: which results in: Let s w;eq denote the saturation where capillary pressure is zero, i.e. P c ðs w;eq Þ ¼ 0. The above equation can be integrated to find the saturation distribution starting from s w ðx ¼ LÞ ¼ s w;eq . The pressure gradients of oil and water at steady state follow from (3) combined with (14): The above corresponds to Darcy's law, where the water flux is constant equal to u T F and the mobilities vary along the core according to the steady state saturation distribution found from (18). System (18) can be solved by separation into a space coordinate integral and a saturation integral: Although the former is trivial, the latter in most cases requires numerical methods. Note that the saturation integral above also has been expressed using normalized saturation S which also facilitates the J ' notation as follows:

Analytical solutions for steady state
We will consider function forms that can solve the above integral analytically. Simple functions with few parameters allow us to illustrate the system behavior and controlling features. We will let both the J-function and the relative permeabilities be linear functions as follows: where J ' < 0 and k * w ; k * o > 0 are all constants. These functions capture that the mobility of a phase increases with its saturation and that the capillary pressure between oil and water decreases with water saturation. Especially, the capillary pressure crosses the saturation axis at a specified value and is given a coefficient J ' to describe the magnitude of capillary forces.
It is convenient to introduce S r , which denotes the reference scaled saturation that is obtained uniformly in the core at steady state if no end effects are present. At steady state, the fractional flow function in the core must be identical to the injection flow fraction. This condition defines the saturation in the core: Applying these functions in (20) gives: It is convenient to introduce the notations: with this notation we observe that the reference scaled saturation in (25) is equivalent to: The integral equation (26) can then be expressed as: The integrand can be reformulated as: which can easily be integrated. The solution is then written in terms of the scaled distance from the outlet Y ¼ y=L as function of normalized saturation and a dimensionless capillary number N: Note that the capillary number N here is defined as a ratio of viscous to capillary forces. From the above equation (31) we see that the scaled position of a saturation from the outlet, YðSÞ, will depend only on the three parameters S eq ; S r ; N. As seen from (31), when N is large, all saturations will be pressed towards the outlet ðY ¼ 0Þ indicating strong advective dominance over the capillary forces. The saturation at Y ¼ 1, denoted S 1 , is of special importance and is defined from (31) by solving: It is important to note that S 1 is not an input, but an output from the solution. Although this equation cannot be solved explicitly, the logarithmic term can be extracted: The average saturation can then be calculated as follows: Note that the logarithmic term (34) has been used to simplify the expression. Next, from (19) and (23) we calculate pressure gradients for each phase as follows: The pressure gradients can be integrated over the core to give the pressure drop of the water phase: and the oil phase:

Comparison with reference behavior (no end effects)
As mentioned, if there are no capillary forces the (uniform) saturation S r is obtained in the core according to equality between the flowing fraction and the injected fraction. Applying this result in (42) and (45) the pressure drops (marked by r to denote reference behavior) become: which are equal to each other. This is expected when capillary forces can be neglected. Both will simply be denoted Δp r . Summarized, the solutions for average saturation and pressure drop can then be expressed as: Further, the difference between the scaled pressure drops can be calculated: When the capillary end effect region is confined within the core we have S 1 ¼ S r which gives the following relations: For this situation, everything within the square brackets is constant and all relations are directly linear with 1=N (or by implication the inverse rate 1=v T ). This proves the intercept method for steady state experiments provided the saturation functions are of the simple form in (23).

Results and discussion
As input parameters we select values representative of water and ndecane, σ ow from Zeppieri et al. (2001) and permeable sandstone (Peksa et al., 2015). J ' ðS eq Þ was reported by Zhou et al. (2017) as À 0.180 for North Sea sandstone and À 0.05 for Liege outcrop chalk, both mixed-wet. Tavassoli et al. (2005) reported a value of À 0.19 for the imbibition curve of strongly water-wet Berea sandstone. An intermediate value of À 0.10 was used. The base rate was 1 pore volume (PV) per day and the base injected water fraction F was 0.1. All parameters are listed in Table 1. From the base parameters, we can calculate characteristic numbers A; B; S r ; S eq ; M; N as listed in Table 2. A comparison of the analytical model with numerical solutions from a commercial software can be found in the Appendix.
Although the model apparently requires specification of all the 14 parameters in Table 1 it is noted that the solution for saturation profile, average saturation and scaled pressure drops only depends on the four parameters M; N; F; S eq while S r and S 1 are intermediate parameters calculated from those four.

Saturation profiles
In this section we show the impact of various parameters on the (normalized) saturation profile SðYÞ following from (31). Using base case parameters the saturation profiles are plotted in Fig. 1 for a high and a low injected fraction; F ¼ 0:9 and F ¼ 0:1, in both cases showing nine profiles corresponding to nine total injection rates varied from 1 PV/d to 100 PV/d with the same factor 10 0:25 � 1:78. The values of the dimensionless number N correspondingly increase from 0.1 to 10.
It is seen that for a given F, at low rates the end effects strongly affect the saturation profiles over the entire core. At higher rates the inlet saturation stabilizes to the value of S r (which is rate independent, see (25)) and even higher rates cause this saturation to cover more of the core. The saturation profiles in all cases converge to SðY ¼ 0Þ ¼ S eq at the outlet. It is also seen that the magnitude of N reflects the extent of end effects. For small values N < 0:1 the terms associated with deviation from S r become large and end effects greatly affect the entire core. For N � 1 (the curves with this value are highlighted) the terms have similar magnitude as the constant terms and the behavior near the inlet is little affected by the end effects. For large values N > 10 the terms representing end effects approach negligible and only a small region near the outlet is affected.
Next we show an example where the oil viscosity is increased to 10 cP, while the other parameters are kept the same as in the previous example. The resulting saturation profiles are depicted in Fig. 2. Notably the higher oil viscosity has led to a lowering of the profiles. This is because a more unfavorable mobility ratio raises the fractional flow function. The saturation giving same fractional flow as the injected fraction must therefore be lower when the oil viscosity is increased. The new values of S r follow from Eq. (27) and are 0.47 for F ¼ 0:9 and 0.011 for F ¼ 0:1.
From (32) the dimensionless number N is seen to be proportional to A þ B. This sum, but not the remaining factor of N, called N 0 , see (27), depends on flow fraction and viscosities. When the two fluids have different mobilities, more viscous force will follow by increasing the fraction of the less viscous (lower mobility) fluid. This increases the magnitude of N and for the same rates of 1-100 PV/d the range of N now corresponds to 0.2 to 20 for the high fraction F ¼ 0:9 and 0.94 to 94 for Table 1 Reference input parameters.  Table 2 Characteristic numbers corresponding to the reference parameters. the low fraction F ¼ 0:1. This high increase compared to 0.1 to 10 for the base case is most significant for the low fractions since the flow is dominated by a greater content of low mobility fluid (oil in this case). Consistent with the description of the previous example, this shift has led all the profiles of the low fraction case to have end effects limited to only a portion of the core while the inlet (and almost half the core for the lowest rate case) obtain saturations equal to S r , not affected by the end effects. For the high injected fraction where the range of N still covers values significantly below 1 the saturation profiles still include cases with all saturations along the core differing greatly from S r . In both cases the value N � 1 seems to reflect a state where a significant part of the core is affected by end effects, while a significant part is little affected ðS � S r Þ. For a given injected fraction, the greatest extent of end effects is seen for the lowest total injection rate of 1 PV/d, corresponding to N ¼ 0:2 for F ¼ 0:9 and N ¼ 0:94 for F ¼ 0:1.

Factors influencing the inlet saturation
The inlet saturation S 1 is indicative of the extent of end effects and is the first saturation to approximate S r , the saturation representative of no end effects. When S 1 � S r the end effect region is limited to within the core and the intercept method becomes valid. It is critical when this value is obtained since it indicates when representative values of the measurement exist in parts of the core.
In Fig. 3 we show the inlet saturation S 1 plotted as function of capillary number N for different mobility ratios M ¼ 0:01; 1; 100 with fractions F ¼ 0:1 and F ¼ 0:6 where the left figure assumes S eq ¼ 0:7 (mixed-wet) and the right assumes S eq ¼ 0:99 (strongly water-wet). For ease of comparison the results are plotted relative to S eq and S r with the fraction S1À Seq SrÀ Seq such that a value of 0 indicates S 1 ¼ S eq (completely capillary controlled) and 1 indicates S 1 ¼ S r (this point is unaffected by end effects). In all cases it is seen that when N passes � 1, S 1 � S r . For S eq ¼ 0:7 the dependence of F and M is not strong and the curves overlap to great extent. S r is obtained for 0:3 < N < 1:5. For the case with S eq ¼ 0:99 there is more spread between the curves and especially mobility ratio matters in addition to N. It should be noted that for the combinations of parameters with high F and low M, a high value of S r is obtained. Particularly, for M ¼ 0:01 we obtain S r ¼ 0:9934 for F ¼ 0:6 and S r ¼ 0:9174 for F ¼ 0:1. These cases that deviate most from the trend with N hence also appear to be cases where end effects matter the least as they have relatively narrow saturation intervals (ΔS ¼ 0:0034 and 0:0726, respectively).

Average saturation and pressure drops
The inlet saturation is necessary input for calculation of average saturation and phase pressure drops according to formulas (48), (49) and (50). Based on these formulas we plot S À S r ; Δpw Δp r ; Δpo Δp r vs 1 N in Fig. 4, Fig. 5 and Fig. 6, respectively and compare with the linear equations (52)-(54) obtained by the assumption of S 1 ¼ S r .
In the three figures it is seen that when 1 N →0 the three parameters follow a linear trend with 1 N in line with the intercept theory. Especially, they overlap with the lines where we have set S 1 ¼ S r . The transition to linear behavior takes place for 1 N � 1:5 to 5 for the different cases, consistent with when S 1 approaches S r as discussed in the previous example.
As seen in Fig. 4, the average saturation with end effects can be higher or lower than the saturation without end effects. If S eq > S r the saturation distribution will be lifted towards higher values than S r , and opposite. Similarly, in Figs. 5 and 6 we note that the pressure drop of a given phase with end effects can both be higher or lower than if end effects were not present. Mainly, end effects will cause a nonuniform saturation and mobility distribution for a given phase. When the end effects shift phase saturations up towards higher phase mobility the Fig. 2. Spatial distributions of scaled saturation S against scaled distance from outlet Y for low ðF ¼ 0:1Þ and high ðF ¼ 0:9Þ injected fractions where the total injection rate has been varied by factors of 10 0:25 from reference (1 PV/d) and upwards (until 100 times the base rate). Higher rate corresponds to higher N with increasing direction indicated by the arrows. The curves with N � 1 are highlighted with bold. Fig. 3. Inlet saturation S 1 (relative to S eq and S r ) plotted vs capillary number N for different mobility ratios (0.01, 1 and 100) and injected fractions (0.1 and 0.6). The left figure shows results for S eq ¼ 0:7 and the right for S eq ¼ 0:99. Fig. 4. Average saturation minus S r plotted against inverse capillary number 1=N for S eq ¼ 0:7 and different mobility ratios using injected fractions F ¼ 0:1 (left) and F ¼ 0:6 (right). Full lines correspond to the analytical solutions, while the dashed lines correspond to the intercept method.

Fig. 5.
Scaled water pressure drop Δp w =Δp r plotted against inverse capillary number 1=N for S eq ¼ 0:7 and different mobility ratios using injected fractions F ¼ 0:1 (left) and F ¼ 0:6 (right). Full lines correspond to the analytical solutions, while the dashed lines correspond to the intercept method. Fig. 6. Scaled oil pressure drop Δp o =Δp r plotted against inverse capillary number 1=N for S eq ¼ 0:7 and different mobility ratios using injected fractions F ¼ 0:1 (left) and F ¼ 0:6 (right). Full lines correspond to the analytical solutions, while the dashed lines correspond to the intercept method. pressure drop will be lower, and when they shift phase saturations down towards lower phase mobility the phase pressure drop will be higher.
At the other extreme, when N→0 the total injection rate is negligible. Setting v T ¼ 0 in (18) gives d x s w ¼ 0, which together with the fixed capillary pressure condition results in SðxÞ ¼ S eq . The scaled pressure drops at this condition naturally correspond to the mobilities at this saturation and defines the stable values that would be approached by the full lines in Figs. 5 and 6 if extended to 1 N →∞. Finally, we note that the difference in the phase pressure drops (for given M; N; F; S eq ) corresponds exactly to the capillary pressure at the inlet, see (51). If S 1 > S eq this corresponds to a negative capillary pressure (and a negative capillary pressure distribution in the core), while if S 1 < S eq a positive capillary pressure is obtained at the outlet and a positive capillary pressure distribution in the core. Whether or not the capillary pressure at the inlet is positive or negative determines which phase has the highest pressure at the inlet. This can have implications for interpreting experimental data as usually one pressure drop is reported from standard experimental designs and one should then be aware of which phase pressure is measured.

Conclusions
In this paper we have derived an explicit analytical solution for coinjection of immiscible phases in a porous medium with capillary end effects at steady state, a setup representative for measuring relative permeabilities. Explicit expressions of saturation profile, phase pressure gradient profiles, average saturation and phase pressure drops were obtained.
� A new capillary number termed N was derived incorporating core and fluid parameters, but even more interesting; also the saturation function parameters, flow fraction and mobility ratio. � The magnitude of the capillary number was very characteristic of end effect behavior. For N < 0:1 end effects were dominant across the entire core; when N � 1 the inlet saturation obtained the saturation representative of no end effect behavior; and when N > 10 most of the core displayed saturations representative of no end effect behavior. � The intercept method was derived analytically. This method states that average saturation and pressure drop (scaled by no end effect pressure drop) can be plotted linearly against the inverse of total injection rate (inverse of capillary number) towards the theoretical values unaffected by capillary end effects. Our method also derives the slope, and the behavior of these parameters in the region where the intercept method (linear behavior) is not valid. The critical value of N where the intercept method could be applied was approximately 1.

Declaration of competing interests
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix. Comparison with a commercial simulator
In this appendix we show how the analytical solutions compare against numerical solutions from the commercial core scale simulation software ProRes Sendra v2018.2.5. The reference parameters in Table 1 were assumed except that six total injection rates were applied; from 1 PV/d and increased by factors of 2 up to 32 PV/d. The software was run with sufficient time for each total injection rate that steady state was achieved. The scaled saturation profiles are shown in Fig. 7 with circles representing the analytical solution and full lines the numerical solution. Near to perfect overlap is obseikrved as expected. Better match could be obtained by tuning the accuracy of the numerical solution. Fig. 7. Comparison of the analytical solution (circle points) with numerical solutions (full lines) generated by a commercial simulator Sendra. Scaled saturation profiles are presented and the reference case input parameters have been used. Six injection rates from 1 PV/d increased by factors of 2 up to 32 PV/d have been applied.