Advances in Water Resources The impact of capillary backpressure on spontaneous counter-current imbibition in porous media

We investigate the impact of capillary backpressure on spontaneous counter-current imbibition. For such displacements in strongly water-wet systems, the non-wetting phase is forced out through the inlet boundary as the wetting phase imbibes into the rock, creating a ﬁnite capillary backpressure. Under the assumption that capillary backpressure depends on the water saturation applied at the inlet boundary of the porous medium, its impact is determined using the continuum modelling approach by varying the imposed inlet saturation in the analytical solution. We present analytical solutions for the one-dimensional incompressible horizontal displacement of a non- wetting phase by a wetting phase in a porous medium. There exists an inlet saturation value above which any change in capillary backpressure has a negligible impact on the solutions. Above this threshold value, imbibition rates and front positions are largely invariant. A method for identifying this inlet saturation is proposed using an analytical procedure and we explore how varying multiphase ﬂow properties af- fects the analytical solutions and this threshold saturation. We show the value of this analytical approach through the analysis of previously published experimental data.


Introduction
Spontaneous imbibition refers to the process by which a wetting phase displaces a non-wetting phase in a porous medium under the action of capillary forces ( Morrow and Mason, 2001 ). Spontaneous imbibition is of particular importance in naturally fractured reservoirs. Carbonate reservoirs, most of which are naturally fractured, contain more than 60% of conventional crude oil reserves ( Beydoun, 1998;Montaron, 2008 ). Hydrocarbons contained within the matrix-blocks of such reservoirs are produced under the action of capillary forces, as water in the fractures imbibes into the matrix-blocks causing the expulsion of the resident oil. Depending upon the boundary conditions, one of three flow mechanisms may occur: counter-current flow, co-current flow or a combination of the two. We define this terminology below.
Consider a core plug that has connate water saturation (the wetting phase) and is otherwise saturated with oil (the nonwetting phase). Suppose that the rock sample is water-wet and that one face of the core plug is exposed to brine while all other faces of the plug are blocked, that is to say that all other boundaries in the system are closed. As the brine spontaneously imbibes * Corresponding author. through the open face upon contact, oil is produced at the water inlet at an exact and opposite flow rate to that of the brine. The total flux is therefore zero. This flow mechanism is counter-current flow ( Nooruddin and Blunt, 2016 ). Now consider the same core plug described above, only this time with two opposite open faces, where one open face, the inlet, is exposed to brine as before, while the second open face, the outlet, is exposed to oil. Suppose further that we introduce a semipermeable membrane, permeable only to the wetting water phase, at the inlet preventing the backflow of oil ( Chen et al., 1992 ). As the water imbibes into the core plug upon contact, the oil flows in a unidirectional fashion to the water. This is usually termed cocurrent flow ( Nooruddin and Blunt, 2016 ).
Removing the porous plate results in oil production at both the inlet and the outlet simultaneously as we observe co-current flow with some degree of counter-current flow. This process is termed free spontaneous imbibition ( Dong et al., 1998 ). The free spontaneous imbibition displacement process can be simulated in the laboratory by exposing one face of the core plug to brine while simultaneously exposing the opposite face to oil, in an otherwise closed system ( Haugen et al., 2014;Mason and Morrow, 2013;Ruth et al., 2015 ).
Several mathematical models describing spontaneous imbibition have been presented in the literature ( Chen, 1988 Fokas and Yortsos, 1982;Kashchiev and Firoozabadi, 2003 ). However, these solutions assumed specific forms for the relative permeabilities and capillary pressure governing the displacement ( Bjørnarå and Mathias, 2013 ). McWhorter and Sunada (1990) presented a semi-anayltical approach for arbitrary saturationdependent flow functions. Schmid et al. (2011) showed that this is a general solution for spontaneous imbibition and used it to derive scaling groups to quantify the imbibition rate ( Schmid and Geiger, 2012 ). The solution, however, is dependent upon the imposed boundary conditions, in particular, the water saturation at the inlet. How one specifies this boundary condition is a question which, until now, has not been properly addressed. A related problem is the so-called capillary backpressure, or the capillary pressure at the inlet of the core. For a core in contact with a free body of water, it has been assumed that this capillary pressure is zero ( McWhorter and Sunada, 1990;Schmid et al., 2011 ). However, for a strongly water-wet rock, this implies that the non-wetting phase is disconnected, and hence it cannot flow: instead there must be a finite capillary pressure, greater than or equal to the entry pressure, to allow the non-wetting phase to escape; this is visually evident in experiments, where the non-wetting phase appears at the face of the core in the form of small bubbles of droplets, with a clearly finite curvature and hence capillary pressure ( Mason and Morrow, 2013 ). In their original mathematical development, McWhorter and Sunada (1990) stated that the inlet saturation is dependent upon the water flux at the inlet. They argued that the appropriate way to handle this from a physical perspective is to apply some flux at the inlet and then implicitly determine the corresponding inlet saturation. In their review, Mason and Morrow (2013) suggested that the saturation applied at the open face of the medium is a serious problem since each saturation value yields a unique solution. In addition, they hypothesized that the correct saturation to apply at the inlet is that yielding the fastest advance of the water as the corresponding solution overwrites all others. Arabjamaloei and Shadizadeh (2010) also proposed that the inlet saturation giving the largest imbibition rate should be used. The truly correct inlet saturation to be applied is still unclear, and hence so is the impact of capillary backpressure since, at a macroscopic level, it is a assumed to be mainly a function of this inlet saturation.
The solution presented by McWhorter and Sunada (1990) is valid only for counter-current and co-current imbibition. For free spontaneous imbibition, a simplified model was presented by Haugen et al. (2014) , where capillary backpressure is incorporated into the theory, although the solution is less general as it assumes piston-like displacements and hence the impact of saturation distributions inside the medium is ignored. Nooruddin and Blunt (2016) extended the analytical solution presented by McWhorter and Sunada (1990) to free spontaneous imbibition as well as other general flow mechanisms where the degree of counter-current and co-current flow within a given system varies with time, although the capillary backpressure was ignored.
In some cases, the non-wetting fluid may be produced as droplets that take some time to form and detach, making the capillary backpressure a function of time, as has been observed experimentally ( Unsal et al., 2009 ). However, most of these experiments were conducted for simple pore systems consisting of a few capillary tubes. In actual rock systems that are sufficiently large, with a few hundred thousand interconnected pores, experiments have indicated production of non-wetting fluid across the entire boundary face upon immersion ( Li et al.,20 09;20 03;20 06 ). This suggests that continuum modelling with a constant backpressure is applicable. Furthermore, these experimental results show a square-root-of-time dependence of production for counter-current imbibition as predicted in the self-similar solution of McWhorter and Sunada (1990) . Therefore, the dynamical behaviour of the system can be modelled using the multiphase Darcy law, the continuum capillary pressure, and conservation of volume, as used commonly to describe the average behaviour of rocks at large scales. The analytical solution assumes that the capillary backpressure is not a function of time and capillary backpressure is uniformly distributed across the inlet boundary.
In this paper we will present analytical solutions for different values of the inlet saturation in counter-current imbibition, corresponding to different amounts of capillary backpressure. For a strongly water-wet rock with a single open boundary, this will be the lowest capillary entry pressure necessary for the non-wetting phase to flow. In weakly water-wet or mixed-wet rocks, the nonwetting phase has a finite relative permeability at zero capillary pressure, and so the capillary backpressure is also zero. However, it is possible to conceive of an experiment where a thin porous plate is attached to the inlet. Imagine that the porous plate is water-wet: in this case the capillary backpressure will be the entry pressure for the non-wetting phase to flow through this plate. We could consider performing experiments using porous plates with different entry pressures: this reproduces the mathematical exercise we perform below.

Governing equations
The treatment follows McWhorter and Sunada (1990) and Schmid et al. (2011) and is repeated here briefly for completeness: the new aspect of the analysis is to explore the sensitivity of the solutions to the inlet saturation (or backpressure). Let q w ( x, t ) denote the water flux at some distance from the medium inlet, x , at some time, t. S w is the water saturation and φ the porosity. It follows from the principle of mass conservation for incompressible flow that Using Darcy's law and assuming counter-current flow with a zero total velocity we find:  where we define the non-linear diffusion term, D ( S w ), as where k is the absolute permeability, μp is the mobility of phase p, k rp is its relative permeability, μ p is the phase viscosity and p c ( S w ) denotes the capillary pressure, defined as the difference between oil pressure and water pressure: where p o and p w denote the oil and water pressures, respectively. We assume that p c is a function of S w only. Now, combining material balance ( Eq. 1 ), with Eq. 2 , we obtain We find a solution as a function only of the variable ( Fu čík et al., 2007 ): Introducing this transformation variable, λ, yields the boundary where we denote q w (x = 0 , t ) = q x 0 (t) . The parameter C is termed the imbibition rate parameter and may be interpreted as a measure of the ability of the porous medium to imbibe the wetting fluid ( McWhorter and Sunada, 1990 ).
We define the fractional flow, F w (S w (x, t )) = F w (x, t ) , accounting for viscous and capillary forces, as Note that F w ( S w ) should not be confused with the commonly used f w ( S w ) which accounts for viscous forces only. We define functions G and T as respectively. The expressions for displacement measured from the inlet of the porous medium, x , the imbibition rate parameter, C , the capillary-driven water fractional flow, F w ( S w ), and its first total derivative with respect to water saturation, F w (S w ) , are given by The derivation of these expressions are provided in Appendices A and C . Capillary backpressure , which is our topic of interest, is defined as the capillary pressure at the inlet of the porous medium, that is In other words, the capillary backpressure is a function of water saturation at the inlet of the porous medium, i.e. p c (x = 0 , t ) = p c (S x 0 ) , as long as the continuum modelling assumption is valid.
It is then clear from Eq. 13 that F w ( S w ) is dependent on the selection of S x 0 , and hence dependent on capillary backpressure. The capillary backpressure is used to define a boundary for the problem.

Methodology
We analyse the effects of varying the imposed water saturation at the inlet boundary, S x 0 . We vary the normalized inlet water saturation while taking a fixed set of relative permeability, k rp ( S w ), and capillary pressure, p c ( S w ), values.
In our computations, we make use of the governing equations given in Section 2 , and the iterative algorithm outlined by McWhorter and Sunada (1990) ; Nooruddin and Blunt (2016) ; Schmid et al. (2011) namely Eqs. 11, 12, 13 and 14 , respectively. This algorithm is described in Appendix D .
Let x f denote the distance advanced by the water front. Mathematically speaking, this is the limit of the distance from the inlet, x ( S w , t ), as the water saturation, S w , approaches the connate water saturation, S wc . Hence x f is defined as We introduce the normalized solution (C D , x f D ) , where the dimensionless imbibition rate and advance of the water front are defined as respectively, for each saturation, S x 0 . For any given inlet satura- We define a threshold water saturation , S * m , as the smallest water saturation at the inlet above which the imbibition rate parameter, C D , and the water front displacement, x f D are effectively invariant. For all practical purposes, we require that C D and x f D are sufficiently close to their respective maximum values. More precisely, we will say that S * x −1 respectively, f −1 denotes the inverse of some function f , = 1 − τ and τ is some arbitrarily chosen tolerance. In this study we set τ to equal 10 −2 . C D and x f D are monotonically increasing in S x 0 . This can be seen immediately from Eqs. 11 and 12 by recognising that they are both The effect of increasing the capillary exponent, n c , is a slower imbibition rate. Note that in (e) the curve corresponding to n c = 1 has a shape differing to the others. This is due to the constant gradient of the capillary pressure curve that can be seen in (a) for is monotonically increasing in S x 0 and it follows that C D and x f D are also monotonically increasing in S x 0 . Also note that in general S * m 1 and S * m 2 are not necessarily equal, that is to say the inlet saturation producing the largest imbibition rate parameter, C D , will not necessarily be the inlet saturation producing the largest water front displacement, x f D . However, the largest amongst the two values will always satisfy Eqs. 20 and 21 . This follows from the definition of S * m , and the fact that C D and x f D are monotonically increasing in S x 0 .

Illustrative example
We assume the following functional forms for relative permeabilities for oil and water phases respectively: as well as the following expression for the capillary pressure, Differentiating with respect to S w gives dp c ( S w ) The other parameters used are listed in Table 1 . The water saturation used in this and the next sections represent a normalized saturation S n :  For illustrative purposes in Table 1 S wc = S or = 0 , and hence Eq. 26 reduces to S n = S w ; in reality for a water-wet rock, typical values of S wc and S or will be around 0.3 ( Blunt, 2017 ). Using these fixed synthetic data; the assumed relationships between relative permeability, capillary pressure and water saturation; and the results presented in Eqs 11, 12 , and 13 ; we vary S x 0 incrementally by 1 100 over the interval [0.01, 1] and compute so- The data is then normalized as previously described, so that we obtain the dimensionless solutions (C D , x f D ) , which are then plotted in Fig. 1 . There exists some threshold inlet water saturation, S * m < 1 , such that C D and x f D are little changed for all S x 0 S * m . From the plot, this value for S * m appears to be approximately near a normalized value of 0.8. This suggests that varying water saturation at the inlet and hence capillary backpressure has minimal impact on C D and x f D at a normalized water saturation of S x 0 = S * m 0 . 8 and greater, where S * m denotes the threshold water saturation at the inlet for this particular choice of parameters.
Fixing the water saturation at the inlet boundary, S x 0 , to equal the maximum mobile water saturation, 1 − S or , will give us the largest imbibition rate parameter, C D , and water front displacement, x f D .
As mentioned previously, Mason and Morrow (2013) stated that varying water saturation at the open face while fixing all other parameters yields a unique solution for each imposed saturation, S x 0 . While this is mathematically correct, from a practical perspective, changes in the inlet saturation have little impact on the solution within a certain range. While there exists a unique solution We can obtain a very close approximation to (C D , x f D ) for a range of inlet saturations and capillary backpressures by fixing S x 0 = 1 − S or . In this example, the choice of which saturation to apply at the open face is not a serious problem at all, assuming that the inlet saturation is above S * m . Fig. 2 shows how the water saturation propagates into the porous medium, when the water saturation at the inlet, S x 0 , is varied. We notice that while the speed of the water front increases when S x 0 is increased from 0.5 to S * m 0 . 8 , increasing S x 0 any further has little effect. The plot corresponding to S * m 0 . 8 overlays that corresponding to S x 0 = 1 . Fig. 2 illustrates how the fractional flow, F w ( S w ), changes with water saturation, S w , and how this behaves as the inlet saturation, S x 0 , is varied. Once again, while an increase in S x 0 from 0.5 to S * m 0 . 8 causes the fractional flow, F w ( S w ), to shift towards the right along the S w axis, this does not occur when we increase S x 0 beyond the threshold inlet saturation, S * m 0 . 8 . We see the respective plots for F w ( S w ) corresponding to S * m 0 . 8 and S x 0 = 1 tracing the same trajectory. Setting S x 0 > S * m results in essentially no change in F w ( S w ). The reason for this can be seen from Eq. 13 , since F w ( S w ) is dependent on T ( S w ), which is in turn dependent on G ( S w ).
The threshold value of 0.8 is not universal as it varies dependent on the relative permeabilities and capillary pressure as we demonstrate in Section 5 . Furthermore, this is a normalized saturation: the actual value will always lie below 1 − S or , where S or is the residual oil saturation.

Sensitivity analysis
We now perform a sensitivity analysis on our results through varying four parameters: end-point mobility ratio, M , capillary pressure, p c , oil relative permeability, k ro , and water relative permeability, k rw . The analysis is performed by varying one parameter at a time while the other parameters are held fixed. The functional forms for k ro , k rw , and p c in Eqs. 22, 23 , and 24 are used, respectively. When they do not vary, the parameters listed in Table 1 are used, except for n c , where a value of 2 is used throughout this section. (In Section 4 a value of n c = 3 is used).

Sensitivity to end-point mobility ratio
Starting with the definition of end-point mobility ratio, M , we have In our illustrative example, our relative permeability end-points are k ro (S wc ) = 1 and k rw (1 − S or ) = 1 2 as in Eqs. 22 and 23 , respectively. Hence, we have and so it is clear that if we vary M , while fixing the water viscosity, μ w , we are effectively varying the oil viscosity μ o .
Here, we change mobility ratio by changing viscosity ratios. The impact of altering the end-point relative permeability values should be similar, provided that these end-points are scalar multiples of their respective relative permeability functions that can be factored out without altering the shape of the curves. This condition is valid for our functional forms used to generate relative permeability as defined in Eqs. 22 and 23 .
We vary the end-point mobility ratio, M , over the set of values {0.05, 0.5, 5, 50}. Results are summarized in Figs. 3 and 4 . The saturation profiles are normalized as illustrated previously in Eq. 26 .
A smaller end-point mobility ratio, M , corresponds to a faster imbibition rate. This is different from the Buckley-Leverett theory for water-flooding where the rate of fluid advance is fixed by the injection rate, although the shape of the water saturation does depend on mobility ratio.
For each of the values we have taken for M , there exists some threshold inlet saturation, S * m < 1 , such that the integral,  the curve. We observe that an increasing end-point mobility ratio, M , corresponds to a decreasing threshold inlet saturation, S * m . By definition, G ( S w ) is dependent on the diffusion parameter, D ( S w ). As such, we perform a similar analysis on D ( S w ). In Fig. 3 c we observe that the behaviour of D ( S w ) is analogous to that of G ( S w ) for any given end-point mobility ratio, M i . The reason for this follows immediately from the definition of G ( S w ) and is to be expected. Fig. 4 show how solutions for C D and x f D behave with varying end-point mobility ratio, M . Note that for each M there exists a threshold inlet saturation, S * m , such that (C D , x f D ) remains almost constant for all S x 0 S * m . This is completely consistent with our findings for G ( S w ) where we tested for sensitivity to M . Once again, we observe that an increasing M corresponds to a decreasing threshold inlet saturation, S * m .

Sensitivity to capillary pressure
The sensitivity of capillary pressure is demonstrated by varying the capillary pressure exponent n c in Eq. 24 over the set of values {1, 2, 5, 10} as shown in Figs. 5 and 6 .
As seen previously in our analysis for end-point mobility ratio, M , there is some S * m such that G ( S w ) is negligible for all S w S * m , for each value of n c . D ( S w ) shows behaviour analogous to that of G ( S w ) when varying n c . An increasing n c corresponds to an increasing water fractional flow F w , at a given S w . The saturation profile shows that the effect of increasing the capillary exponent, n c , is a slower advance of the water. Note that in Fig. 5 e the curve corresponding to n c = 1 has a shape differing to the others. This is mainly due to the constant gradient of the capillary pressure curve, as indicated by Eq. 25 for n c = 1 .
We now demonstrate the effects on C D of varying n c as shown in Fig. 6 a. For each n c there exists an S * m < 1 such that C D remains essentially unchanged for all S x 0 S * m . In Fig. 6 b we see that for each n c , there exists an S * m such that x f D is practically constant for all S x 0 S * m . We see that increasing n c corresponds to decreasing S * m for both C D and x f D , consistent with our findings for G ( S w ) and D ( S w ).

Sensitivity to oil relative permeability
The sensitivity to k ro is demonstrated by varying the oil exponent n o in Eq. 22 over the set of values {1, 2, 5, 10}. The results are shown in Figs. 7 and 8 . Fig. 7 a shows the fixed capillary pressure curve used for this analysis. By varying the oil exponent, n o , we obtain the set of oil relative permeability curves presented in Fig. 7 b. As seen previously in our analyses for end-point mobility ratio, M , and capillary pressure exponent, n c , Fig. 7 f shows that there is some S * m such that G ( S w ) is negligible for all S w S * m , for each value n o takes. In Fig. 7 c, D ( S w ) shows behaviour analogous to that of G ( S w ) when varying n o . Fig. 7 d illustrates that an increasing n o corresponds to an increasing water fractional flow, F w , for a given S w . We also see in Fig. 7 e that the effect of increasing the oil exponent, n o , is a slower rate of advance of the water.
The effects on C D of varying n o are shown in Fig. 8 a. For each n o there exists an S * m < 1 such that C D remains essentially unchanged for all S x 0 S * m . We see that for each n o , there exists an S * m such that x f D is practically constant for all S x 0 S * m , as shown in Fig. 8 b. We also see that increasing n o corresponds to decreasing S * m for both C D and x f D once again, consistent with our findings for G ( S w ) and D ( S w ). In comparison with other parameters, it is apparent that S x 0 is most sensitive to k ro .

Sensitivity to water relative permeability
We closely follow the previous procedures to test the sensitivity of k rw by varying the water exponent n w in Eq. 23 over the range of values {1, 2, 5, 10} as depicted in Fig. 9 b.
As seen in previous analyses, Fig. 9 f shows that there is some S * m such that G ( S w ) is negligible for all S w S * m , for each value n w takes. Fig. 9 d illustrates that an increasing n w corresponds to a decreasing water fractional flow, F w , for a given S w . One remarkable feature of varying k rw is demonstrated in Fig. 9 e, in which a sharp frontal displacement behaviour is obtained by increasing n w . This is because we have a very low water relative permeability over much of the saturation range, which limits the rate of imbibition.
We now study the effect of k rw on S x 0 , see Fig. 10 . For each n w there exists an S * m < 1 such that the imbibition rate, C D , remains essentially unchanged for all S x 0 S * m . Similarly for each n w , there exists an S * m such that the position of the water front, x f D , is practically constant for all S x 0 S * m . We see that increasing n w corresponds to decreasing S * m for both C D and x f D once more, consistent with our findings for G ( S w ) and D ( S w ), although the variation here is small compared with that for the other parameters studied previously. shows behaviour analogous to that of G ( S w ) when varying n w . (d) Increasing n w corresponds to a decreasing water fractional flow, F w , for a given S w . (e) Increasing the water exponent, n w , imbibition is slower. (f) There is some S * m such that G ( S w ) is negligible for all S w S * m , for each value n w takes.

Comparison with experimental data
In this section, we compare our methodology with experimental data published by Li et al. (2006) . The rock and fluid properties for the sample used in this analysis are given in Table 2 .
Since the data does not contain independent measurements of capillary pressure and relative permeabilities, Li et al. (2006) used a numerical modelling approach to produce this information by matching experimental data. The functional forms of these functions are given by ( Li et al., 2006 ); where p x 0 c is the capillary pressure at the inlet boundary, p x f c is the capillary pressure at the front, and S x f is the water saturation at the front. Similarly for Eqs. 30 and 31 , the superscripts x 0 and x f are used to indicate the value at the inlet boundary and at the front, respectively. Although not explicitly mentioned, these parameters solely depend on the saturation values at these positions. The parameters in Eqs. 29, 30 and 31 , are therefore alternative representations of saturation end-points.
An analytical approach will be used to produce the multiphase flow functions that match the experimental measurements. Using the functional forms in Eqs. 29, 30 and 31 , capillary pressure and relative permeability curves are generated. An analytical solution is used to compute x ( S w , t ), C, F w ( S w ), and F w (S w ) using Similarly for each n w , there exists an S * m such that the location of the water front, x f D , is practically constant for all S x 0 S * m . We see that increasing n w corresponds to decreasing S * m for both C D and x f D once more, consistent with our findings for G ( S w ) and D ( S w ). 3.14 × 10 3 [Pa] Table 3 Parameters used in Eqs. 29, 30 and 30 to generate capillary pressure, relative permeability to oil and relative permeability to water curves, respectively, for Sample # H8O. (Modified from Li et al. (2006) ).
Parameter Li et al., (2006) where S w is the average water saturation in the system, ˆ x f D is the normalized water saturation at the front as defined here in Eq. 33 , not to be confused with the dimensionless x f D defined in Eq. 19 , and p o L is the oil pressure at the dead-end boundary. The derivation of Eq. 34 is provided in Appendix E . The set of curves published by Li et al. (2006) will also be used in our analytical procedure and compared with the measured data.
Results are displayed in Fig. 11 , in which experimental data are compared to model outputs with two sets of capillary pressure and relative permeability curves. The first set of solutions uses the properties presented in Li et al. (2006) . Note that the match to the experiments is poor. Li et al. (2006) used a numerical solution that appeared to reproduce the experiments closely; however, we suggest that their numerical solutions suffered from significant discretization errors. As shown in previous work, a very refined grid is needed to capture the imbibition rate accurately ( Nooruddin and Blunt, 2016 ). Instead, a very different set of parameters -see Table 3 and Figs. 11 and 12 -are needed to match the experiments accurately. This illustrates the value of the analytical solution, as it removes numerical errors from the analysis.   Table 3 .
The multiphase flow parameters produced in this study ( Fig. 12 b) indicates strongly water-wet behaviour. The residual oil saturation is 0.4. Furthermore, the end-point relative permeability values are chosen to be consistent with this wetting condition. The capillary pressure curve ( Fig. 12 a) demonstrates a value of capillary backpressure of 3.15 × 10 3 Pa, which is approximately 6% of the capillary pressure at the front. Other important information is also shown in Fig. 13 .
We now estimate the value of water salutation at the inlet boundary ( S x 0 ) using the methodology that was previously described. For this, we compute x f D and C D as a function of S x 0 , as Fig. 13. Computations of (a) capillary diffusion, (b) fractional flow and (c) invariant water saturation profiles using the capillary pressure and relative permeability curves shown in Fig. 12 .   Fig. 14. Computations of x f D and C D as a function of S x 0 . It can be seen that above a value of 0 . 53 − 0 . 54 , the imbibition rate and normalized front location are effectively invariant, suggesting a value of the water saturation at the inlet face boundary.
shown in Fig. 14 . The plot indicates that above a value of approximately 0.53 the solutions remain largely invariant, which is very close to the average water saturation in the core sample at the end of the imbibition test of 0.525.

Conclusions
We have presented analytical solutions for counter-current imbibition with capillary backpressure. Imposing a capillary backpressure is equivalent, at Darcy-scale, to specifying the water saturation at the inlet. The backpressure is either the capillary entry pressure for the non-wetting phase that can be specified using a porous plate, or zero for a rock that is not strongly water-wet with a free boundary, where the non-wetting phase has a finite mobility at the inlet.
We show that for example cases representing a strongly waterwet system that there is a threshold inlet water saturation above which the imbibition rate and the advance of the wetting front are largely invariant. As a consequence, the imbibition behaviour may be independent of the capillary backpressure, assuming that the inlet saturation is close to its maximum value.

Acknowledgments
The second author would like to thank Saudi Aramco for their generous financial support towards his PhD degree at Imperial College London.
Appendix A. Fractional Flow, F w ( S w ), Imbibition Rate Parameter, C , and Distance From the Inlet, x We define the fractional flow, F w (S w (x, t )) = F w (x, t ) , accounting for viscous and capillary forces, as Rearranging Eq. 35 to make q w ( x, t ) the subject and then substituting into Eq. 1 gives Note that q x 0 ( t ) is taken outside of the derivative since it is not a function of x , also note the total derivative dF w (S w ) dS w since F w ( S w ) is only a function of S w . We introduce the transformation variable Substituting q x 0 ( t ) from Eq. 7 and the transformation variable λ, we can rewrite Eq. 37 as where F w (S w ) denotes the total derivative of F w ( S w ) with respect to S w . Hence, by definition of λ, From Eq. 42 λ = λ(S w (x, t )) is only a function of S w , since S w is the only unknown variable. Note also that for capillary driven flow, the displacement of the water front scales as a function of √ t , since t . In the analogous Buckley-Leverett solution for viscous driven flow, the advance of the water front scales as a function of time as Finally, we make a remark on the behaviour of fractional flow at boundaries, namely at the inlet and at the water front. At the inlet, Eq. 44 is a general solution which can be applied to either counter-current or co-current flow. For counter-current flow we substitute q t (t) = q w (x, t ) + q o (x, t ) = 0 into Eq. 44 and the first term on the right-hand side of the equation will vanish. Eq. 44 then becomes Note the total derivative dλ dS w since it follows from Eq. 42 that λ = λ (S w (x, t )) .
Differentiating Eq. 42 with respect to S w gives where F w (S w ) denotes the second total derivative of F w ( S w ) with respect to S w . Combining Eq. 4 8 and 4 9 gives This is a second-order ordinary differential equation with the following boundary conditions for the displacement of the nonwetting oil phase by the imbibing wetting water phase: The physical meaning of Eqs. 51 and 53 is that while the fractional flow at the inlet is equal to 1, the applied water saturation at the inlet, S x 0 , will not propagate into the porous medium, hence its derivative with respect to water saturation is zero. Looking further into the porous medium beyond the inlet, the water saturation, S w , decreases with distance, x , until connate water saturation, S wc , is reached. At this point the fractional flow is zero, giving Eq. 52 . Integrating Eq. 50 twice with respect to S w gives for some constants K 1 and K 2 . For the sake of brevity, we let The derivation for this remark can be found in Appendix B . By making use of the boundary conditions, Eqs. 51, 52 and 53 , it follows from Eq. 54 that The full derivation of these expressions for C, F w ( S w ) and F w (S w ) are provided in Appendix C. F w .
Integrating Eq. 50 with respect to S w gives and integrating once again with respect to S w gives for some constants K 1 and K 2 . For the sake of brevity, we let We make use of Eq. 51 . At S w = S x 0 we have F w (S x 0 ) = 1 and T (S x 0 ) = S x 0 S x 0 G (β ) dβ = 0 , substituting these values into Eq. 66 we obtain 1 = − φ 2 C 2 (0) + K 1 S x 0 + K 2 (67) Next we make use of Eq. 52 . At S w = S wc we have F w (S wc ) = 0 and Eq. 66 becomes Subtracting Eq. 69 from Eq. 68 yields Substituting K 1 into Eq. 68 gives K 2 = 1 − K 1 S x 0 (72) and substituting K 1 and K 2 into Eq. 66 gives Now, differentiating Eq. 75 with respect to S w gives We now make use of Eq. 53 . At S w = S x 0 , we have F w (S x 0 ) = 0 and T (S x 0 ) = s x 0 s x 0 G (β ) dβ = 0 . Consequently, Eq. 76 becomes Eq. 77 implies that 1 − φ 2 C 2 T (S wc ) = 0 , therefore, when substituting in C , Eq. 75 becomes Hence, since T (S w ) = −G (S w ) by Eq. 55 , we also have To determine the advance of the water front at a given water saturation and time, it remains to compute the imbibition rate parameter C . This is done using Eq. 12 , namely C = φ