Capillary Hysteresis and Gravity Segregation in Two Phase Flow Through Porous Media

We study the gravity-driven flow of two fluid phases in a one-dimensional homogeneous porous column when history dependence of the pressure difference between the phases (capillary pressure) is taken into account. In the hyperbolic limit, solutions of such systems satisfy the Buckley-Leverett equation with a non-monotone flux function. However, solutions for the hysteretic case do not converge to the classical solutions in the hyperbolic limit in a wide range of situations. In particular, with Riemann data as initial condition, stationary shocks become possible in addition to classical components such as shocks, rarefaction waves, and constant states. We derive an admissibility criterion for the stationary shocks and outline all admissible shocks. Depending on the capillary pressure functions, flux function, and the Riemann data, two cases are identified a priori for which the solution consists of a stationary shock. In the first case, the shock remains at the point where the initial condition is discontinuous. In the second case, the solution is frozen in time in at least one semi-infinite half. The predictions are verified using numerical results.


Introduction
In this paper we investigate gravity driven flow of two fluid phases through a homogeneous one-dimensional porous column. We are concerned with the special case in which the length of the column is large and no injection of fluid is present (i.e. the total flow is zero). The corresponding mathematical model uses the well-known Buckley-Leverett equation, see [9]. In dimensionless form it reads ∂ t S + ∂ x h(S) = 0. (1.1) Here S ∈ [0, 1] is the saturation of the wetting phase and h : [0, 1] → [0, ∞) the fractional flow function of which a typical sketch is shown in Figure 1 (left). The space coordinate is x and t denotes time. We solve (1.1) for t > 0 and x ∈ R, where we prescribe at t = 0 the Riemann condition S(x, 0) = S T for x < 0, S B for x > 0, where 0 < S B < S T < 1 are constants. (1.2) Solutions of (1.1)-(1.2) are generally non-unique [14,Chapter 1]. To find the physically relevant solution, we replace (1.1) by the capillary Buckley-Leverett equation (1.3) and study the limit δ → 0. In (1.3), p is the pressure difference between the fluid phases and δ is the dimensionless capillary number which scales inversely with the length of the domain (reference length). Hence for large domains, (1.1) is approximated by (1.3). A detailed derivation is given in Section 2.
To solve (1.3) a relation between S and p is assumed. In the standard equilibrium approach, one uses p = p c (S), (1.4) where p c : (0, 1) → [0, ∞) is a capillary pressure function [9]. For δ > 0, let (S δ , p δ ) denote the solution of (1.2)-(1.4). The hyperbolic limit as δ → 0 yields the well known Buckley-Leverett solution of (1.1)-(1.2), comprising of constant states separated by shocks and rarefaction waves [14,15,22]. In particular, shocks can be seen as the limit of smooth, monotone travelling waves which become steeper as δ → 0 [14,15,31]. It is also known that the hyperbolic limit solution is independent of the actual shape of the capillary pressure p c as long as the approximating equation (1.3) is of convection-diffusion type [14,Chapter 3].
However, there are other more realistic capillary pressure expressions for which the limiting (δ → 0) solution does inherit some properties of the vanishing capillarity. This was studied in detail in [19,30,32,34] for the case where (1.4) is replaced by the nonequilibrium expression where τ > 0 is a relaxation parameter called dynamic capillarity coefficient which attributes to saturation overshoot [17,20,33]. In this paper we investigate the effect of hysteresis in the capillary pressure. Since p c (·) is a single valued function of saturation only, it does not contain any information about the history or directionality of the process. In particular, it does not distinguish between imbibition and drainage. No hysteretic effects are present in (1.4). However, hysteresis is known to occur in multi-phase porous media flow. This was first observed by Haines in 1930 and has been verified since by numerous experiments, some notable examples being [21,35]. An overview of different hysteresis models from the mathematical, modelling and physical perspective can be found in [8,16,28]. In our approach we replace (1.4) by the following hysteresis description. Let denote the imbibition and drainage capillary pressure functions, typical examples being shown in Figure 1 (right), and let be the hysteresis region in the (S, p)-plane, as indicated in Figure 1 (right). We restrict ourselves to the well-known play-type hysteresis model, first proposed in [6], which relates S and p via the relation p ∈ H, and p = The play-type model assumes that switching from drainage to imbibition, or vice versa, occurs only along vertical scanning curves. For the rest of the study the state when ∂ t S = 0 and consequently p ∈ [p c (S)] will be referred to as the undetermined state of the system, whereas ∂ t S > 0 (consequently p = p (i) c (S)) and ∂ t S < 0 (p = p (d) c (S)) will refer to the imbibition and drainage states, respectively. Other models, such as the Lenhard-Parker model [23] and the extended play-type hysteresis model [8,18] assume a more complex relation between S and p when p c (S). We comment on the consequences of such models in Remark 4.2.
The purpose of this paper is to construct solutions of the Riemann problem (1.1)-(1.2) that arise as the hyperbolic limit (δ → 0) of (1.2)-(1.3) and (1.5). We demonstrate that the occurrence of hysteresis in the vanishing capillary term, i.e. using (1.5) instead of (1.4), gives solutions that are significantly different specially when S B is close to 0 and S T is close to 1. In particular, a stationary discontinuity occurs at the location where the initial condition is discontinuous. This is illustrated in Figure 2. The magnitude of the jump depends on the difference p  curves are close to each other, whereas in the (right) plot the curves are as shown in Figure 1 (right). Exact details are given later. was mentioned briefly in Section 3.5 of [19]. A stationary discontinuity in saturation for the classical problem (with no hysteresis) was also predicted in [2] at the interface between two semi-infinite halves having different h and p c functions. However, the directionality imposed by hysteresis demands an extension of their results. The saturation discontinuity has major practical importance as the saturation distribution can change considerably if hysteresis is present, as is evident from Figure 2. Remark 1.1 (Vanishing capillarity method). In the mathematical literature, the method of finding solutions of equations such as (1.1) by passing to the limit δ → 0 in (1.3) is called the vanishing viscosity method, and was studied in classical references such as [13,22]. In these papers, the notion of entropy was introduced, calling the vanishing viscosity solutions 'entropy' solutions. In the context of our application, we call the approach vanishing capillarity method, the approximate solutions (S δ , p δ ) for δ > 0 the capillarity solutions, and the δ → 0 limit solution the vanishing capillarity solution.
Hysteresis models have been analysed exhaustively, particularly in one spatial dimension [17]. Existence and uniqueness results for the regularised play-type model are given in [10,12,18,27] and [7] respectively. A horizontal redistribution study using similarity solutions was performed in [8]. Travelling wave analysis for hysteresis was conducted for flow of water through soil in [5,20,33]. For the two-phase case, travelling waves were studied in [19] for monotone flux functions and vanishing capillarity solutions were derived that differ significantly from the classical ones. Non-classical solutions for non-monotone flux functions like h are investigated in [29], although hysteresis is not included. In [11], examples of how hysteresis influences hyperbolic solutions in mechanics are found. The role of relative permeability hysteresis, which is not addressed in this study, in determining the hyperbolic solutions is examined in [1,4,25,26]. However, the relation between p and S is assumed to be a linear one in these articles. In the current study, we investigate the vanishing capillarity solutions for the capillary hysteresis models in the non-monotone flux case and show that non-classical behaviour such as stationary shocks may occur. Stationary discontinuities have been studied in [2] for heterogeneous media without hysteresis and for redistribution problems (no gravity) in [8,24] with hysteresis. The occurrence of a stationary discontinuity due to hysteresis in gravity driven flows is to our knowledge a novel observation.
We structure the paper as follows: In Section 2 the assumptions are stated and the model is derived. Using travelling wave analysis, all admissible shocks, including stationary shocks, are derived in Section 3 for the standard and the hysteresis model. Then, in Section 4, the admissible shocks are used to construct the vanishing capillarity solutions. Two cases are identified when the solutions for hysteresis deviate from the classical ones. This is determined a priori from h, p c , S B and S T . In the first case, the solution has a stationary discontinuity while the rest of the solution retains the structure of the classical solutions (Section 4.2.1). The second case has no classical counterparts and the solution is frozen in one of the semi-infinite halves (Section 4.2.2). Finally, in Section 5, we solve (1.2)-(1.3) and (1.5) numerically for small δ > 0 and show that the solution closely resembles our predictions.

Problem Formulation
We set-up the two-phase flow problem in a one-dimensional homogeneous porous domain. The phases are assumed to be incompressible and immiscible. There is no injection at the boundaries, making the flow purely counter-current and gravity driven, having zero total flux of the combined wetting and non-wetting phases. Following [9], we consider for each phase the mass balance equation and the corresponding Darcy Law. This yields For each phase α = w, n (w and n representing the wetting and the non-wetting phases respectively), p α denotes the pressure, k α the relative permeability, µ α the viscosity and ρ α the density. The porosity φ and absolute permeability K are properties of the medium and are constant due to the assumption of homogeneity. Finally, gravity points in the direction of positive x and g is the gravitational constant. For the remainder of the study we assume that the wetting phase is denser than the non-wetting phase, i.e., Adding the equations in (2.1) gives The term inside the brackets [ ] is the total flux of the combined phases. Since no-injection takes place in the column Defining the capillary pressure and mobility ratio as p := p n − p w and M := µ w /µ n respectively, one obtains through rearrangement To non-dimensionalise (2.1), we introduce a characteristic pressure p ref (taken from the capillary pressure curves), a characteristic length H (length of column or typical observation distance) and a characteristic time t ref . We further define the fractional flow function and introduce the dimensionless capillary number , and redefining the dimensional quantities as their dimensionless versions The Brooks-Corey model is commonly used for the relative permeability functions: This gives the shape of h as shown in Figure 1 and properties outlined in (P2). As for p, either (1.4) or (1.5) is used. For the capillary pressures and the fractional flow function, following set of properties is assumed, see Figure 1. They are consistent with experimental observations [3,9]: The assumption (P1) is consistent with the van Genuchten model for capillary pressures.
For h, the Brooks-Corey model (2.5) yields

Admissible shocks
In this section we consider shock solutions of equation (1.1) that originate from smooth solutions of (1.3). A shock is characterized by a constant left state S l , a constant right state S r and a constant speed c. They are denoted by {S l , S r , c} or We call {S l , S r , c} an admissible (i.e. vanishing capillarity) shock if it can be approximated, as δ → 0, by smooth solutions S δ of (1.3) [14,15]. To investigate this, we consider a special class of solutions of (1.3) in the form of travelling waves: Here S : R → [0, 1] is the saturation profile, p : R → [0, ∞] the capillary pressure profile and c ∈ R the wave-speed. We consider profiles that satisfy for the saturation and for the pressure Clearly, if a smooth profile (3.2) exists and satisfies (3.3), then where primes denote differentiation. Integrating this expression yields If S satisfies (3.3), then (3.6) implies that p ′ has a limit for η → ±∞. The boundedness of pressure (3.4) then forces Using this in (3.6) for η = ±∞ results in Substituting c and A in (3.6) yields 3.1 Equilibrium case: p given by (1.4) Then p = p c (S), and we write (3.10) as Hence, the shock {S l , S r , c} is admissible if c is given by (3.9) and Conditions (3.12) are classical (Oleinik admissibility conditions), see [14,15,22]. To see which possible admissible shocks connect to the left or right states of the Riemann problem Figure 3. Thus, the pointsS B andS T satisfy, Applying (3.12) to Riemann problem (1.1)-(1.2), we note that if S B < S 1 , then S B can serve as the right state S r = S B of an admissible shock for left states between S B < S l <S B , see Figure 3. Similar for S l = S T > S 2 andS T < S r < S T .
3.2 Hysteretic case: p given by (1.5) If c(S l , S r ) > 0 or c(S l , S r ) < 0, then any admissible shock in terms of (1.4) is also an admissible shock in terms of (1.5): To see this, fix S l > S r and let (3.12) hold. Define the relation between p(η) and S(η) as Then, separately for c(S l , S r ) > 0 and c(S l , S r ) < 0, S and p are related through the classical relation (1.4) where p c is replaced by p the existence of S(η) and p(η) solving (3.10) with boundary conditions (3.3)-(3.4) follows from Section 3.1. Observe that S ′ (η) < 0 for all η ∈ R as a consequence of (3.11) and (3.12). Hence, recalling from (3.2) that S δ (x, t) = S(η) and p δ (x, t) = p(η), we have It is now straightforward to verify that (S δ , p δ ) satisfies the hysteresis relation (1.5) for both c(S l , S r ) > 0 and c(S l , S r ) < 0. Hence, shocks that are admissible in terms of the equilibrium capillary pressure (1.4) are admissible in terms of the hysteretic capillary pressure (1.5) as well, provided c = 0. Note that the profiles (S, p) depend on the functional form of p c , but not the resulting shock.
Observe that for c(S l , S r ) > 0, the entire approximating wave (S δ , p δ ) is in imbibition state since p δ (x, t) = p (i) c (S δ (x, t)) and ∂ t S δ > 0. Hence, we also refer to the resulting shock as being in imbibition state. Similarly, for c(S l , S r ) < 0, the shock is in drainage state since p δ (x, t) = p (d) c (S δ (x, t)) and ∂ t S δ < 0. The effects of hysteresis are only observed, if the hysteretic states ahead and behind the shock are different. This can only happen for stationary shocks having zero wave speed, i.e. c(S l , S r ) = 0. For such shocks, Rankine-Hugoniot's expression (3.9) requires (3.14) Note that this implies Without a loss of generality we put the stationary shock at x = 0.

Case A: Connection between imbibition and drainage states
Let us first consider the case The case of S l in imbibition and S r in drainage is symmetrical. Here, the imbibition and drainage states are inherited from the approximating solutions as before.
Repeating this reasoning for η > 0, we conclude that the unique solution of (A1)-(A3) is c (S l )) for η < 0, (S r , p r ) = (S r , p (i) c (S r )) for η > 0, (3.23) where the pair {S l , S r } in addition to satisfying (3.14), also must satisfy p l = p      Connecting an imbibition state with another imbibition state demands that (3.24) be satisfied with p l = p (i) c (S l ) and p r = p (i) c (S r ) which has the unique solution S l = S r . The same result holds for connecting a drainage state with another drainage state. Hence, no non-trivial stationary shock (S l = S r ) exists in these cases.

Case B: One state undetermined
It is also possible that one of the states {S l , S r } is undetermined (neither in imbibition nor in drainage). To demonstrate the admissibility for this case, we first assume that the left state is undetermined and the right state is in imbibition. Using ∂ t S = c(S l , S r )S ′ = 0 in (1.5), we have that {S l , S r , 0} is an admissible stationary shock if (B1) Relation (3.14) is satisfied. Following the same arguments as before, we have (S(η), p(η)) = (S r , p r ) = (S r , p (i) c (S r )) for all η > 0. (3.30) What remains is the η < 0 problem which reads Clearly, S = S l and p = p l is a solution provided, Hence, also in this case, condition (3.24) is satisfied. We show below that additionally The construction is also shown in Figure 5. It is straightforward to verify that this is only case when (3.33) is satisfied.
c (S l )) is also shown.
Condition (3.34) serves as the admissibility criterion when the left state is undetermined and the right state is in imbibition. Considering all combinations, including undetermined-drainage, imbibition-undetermined, drainage-undetermined and undeterminedundetermined, the complete list of admissible stationary shocks {S l , S r , 0} with S l = S r is given in Table 1. Table 1: All admissible stationary shocks {S l , S r , 0} with S l = S r . They satisfy the condition h = h(S r ) − h(S l ) = 0, and p = p r − p l = 0 where S l , S r , p l and p r are given in the table. The states on the left and the right of the shocks are also stated. Imbibition state corresponds to p (i) c curve being used to relate S l and p l or S r and p r . Simialrly, drainage state corresponds to the use of p

Vanishing capillarity solutions of the Buckley-Leverett equation
We construct the solution of the Riemann problem (1.1)-(1.2) with the admissibility conditions from Section 3. A solution of (1.1) is composed of constant states separated by shocks and rarefaction waves. We recall that a rarefaction wave is a smooth solution of (1.1) having the form S(x, t) = r(ζ) where ζ = x/t. where r(ζ l ) = S l , and r(ζ r ) = S r or ζ l = dh dS (S l ), and ζ r = dh dS (S r ).

Rarefaction waves are well-defined if
f ′′ (s) does not change sign for all min{S l , S r } ≤ s ≤ max{S l , S r }. We show below how to construct solutions of (1.1) with rarefaction waves satisfying (4.2) along with condition (4.3), and shocks satisfying the admissibility conditions from Section 3. For completeness, we briefly recall the classical Buckley-Leverett construction.

Classical construction
The construction of vanishing capillarity solutions for the classical case is well established, see [14,22]. Recalling the definition ofS B andS T in (3.13), we identify three cases. Since the solutions for these cases will be used in Section 4.2 to construct parts of the vanishing capillarity solution for pairs other than {S B , S T }, we also use S l and S r in the notation: In this case S B < S M since S B < S T . The vanishing capillarity solutions are (4.4a) and with r(·) defined in (4.2), and (4.5b) Case III: S B < S M < S T : In this case we have,

Construction with capillary hysteresis
Vanishing capillarity solution for Case I and II: Observe that in Case I of Section 4.1 (S T < S M ) the shock is in an imbibition state and the rarefaction wave part, if it exists, satisfies ∂ t S > 0. Hence, the vanishing capillarity solution as a whole is in imbibition state. Similarly, for Case II the entire solution is in drainage state. Consequently, based on the discussions in Section 3.2 and the fact that the classical solution does not depend on the specific form of p c , the vanishing capillarity solutions for the capillary hysteresis model in Case I and Case II are identical to the classical ones.
Vanishing capillarity solution for Case III: As in the classical case, we construct a solution which is in imbibition state when x > 0 and drainage state when x < 0. Hence one expects to have a stationary shock at x = 0. We restrict ourselves to the case the solution for h(S T ) ≤ h(S B ) being symmetrical. With respect to the cases A and B in Section 3.2 we divide the discussion in two parts.

Case
The restriction (4.7) implies that S B ≤ S * . The stationary shock at x = 0 in this case must connect S * and S * . To see this, assume the contrary. Recalling Table 1, we then have lim  Table 1. Thus, the stationary shock at x = 0 connects S * and S * .
Vanishing capillarity solution: For x > 0 and t > 0, S(x, t) is given by (4.4) with S l = S * . Restriction (4.7) then implies that S B <Š T . In this case, the stationary shock {S * , S * , 0} is ruled out. Otherwise, since S T < S * , the classical solution for x < 0 and t > 0 will be in imbibition state. However, recalling Table 1, the shock {S * , S * , 0} is not admissible for left state in imbibition. Similarly, other possibilities are eliminated using Table 1, except the following scenario: Vanishing capillarity solution: for x > 0 and t > 0, S(x, t) is given by (4.4) with S l =Š T . The vanishing capillarity solutions for the cases S T > S * and S T ∈ (S M , S * ) are shown in Figure 6. Observe that there is no classical counterpart of (4.10). The solution plotted in Figure 2 (right) for the play-type hysteresis model is of this type.
Remark 4.2 (Generality of the results with respect to other hysteresis models). The admissible shocks presented in Table 1 and the vanishing capillarity solutions presented in (4.8) and (4.10) are consistent with any hysteresis model that satisfies the condition p ∈ H, and p = implies ∂ t S ≤ 0. (4.11) Most of the commonly used hysteresis models, including the Lenhard-Parker model [23] and the extended play-type model [8] belong to this category. The results are consistent with any model satisfying (4.11) since, the models only differ in the description of (S, p) c (S)), i.e. when the hysteretic state is undetermined. As a result, if a stationary shock connects imbibition to drainage, then the set of equations (A1)-(A3) are valid also for a model satisfying (4.11). Thus, the resulting shocks are unaltered. Similarly, (B1)-(B3) (in particular (B3)) are consistent with describing the shock when one of the states is undetermined since in this case ∂ t S = c(S l , S r )S ′ = 0 which is allowed by (4.11). Hence, models satisfying (4.11) are consistent with admissible shocks listed in Table 1, and consequently with the vanishing capillarity solutions derived in Section 4.

Numerical Results
To solve (1.2), (1.5) numerically, (1.5) is usually regularised in the following way: when ∂ t S < 0. (5.1) The relaxation parameter τ > 0 (dynamic capillarity coefficient) was mentioned briefly in Section 1. The approach for solving the system (1.3), (5.1)-(5.2), is based on rearranging (5.1) to express ∂ t S as a function of S and p, and using this, to consider (1.3) as an elliptic equation for the pressure. Details of the numerical method are given in Section 5 of [19], see also the numerical sections of [20,33]. Cell centered finite difference method with a uniform mesh is used for the computation with ∆x and ∆t representing the mesh and the time step sizes respectively. The following choices of parameters are made τ = 0.01, δ = 0.25, H = 100, ∆x = 0.01, ∆t = 0.0001.
These values ensure both that our parabolic solver is converging and the capillarity solutions are sufficiently close to their hyperbolic limit. Since smaller δ also implies that the profiles take longer to develop, we have optimized the value of δ so that a good approximation of the developed profile is obtained in a reasonable time.  Figure 7 (right). In this case, direct computation shows that S * = 0.4121 and S * = 0.5879. We take S T = 0.8 and S B = 0.1 implying that S T > S * . The vanishing capillarity solution expected in this scenario is outlined in (4.8). It consists of shocks between S T -S T and S B -S B , rarefaction waves betweenS T -S * andS B -S * and a stationary shock from S * to S * at x = 0, see Figure 2 (left) and Figure 7 (left). The numerically computed solution is shown in Figure 7 (left) and in the p c -S plane in Figure 7 (right). The numerically obtained capillary solutions are a close match to the vanishing capillarity solutions predicted. The minor differences are due to δ > 0, τ > 0 and numerical errors. Note that the countercurrent flow for x < 0 stems from non-monotonicity of h since the flux at x = 0 is h(S * ) which is greater than the flux at the left boundary. These curves resemble experimentally obtained retention curves (see [21]) and were taken from [19]. With {S B , S T } same as before, one has S T < S * = 0.8759 in this case. Hence, recalling Section 4.2.2, particularly (4.10), the solution is frozen for x < 0. For x > 0, there is a shock connectingŠ T = 0.2 and S B . This behaviour is mimicked by the computed capillarity solutions, see Figure 8. In this case, the total mass is still balanced since the total rate of water infiltration is c(S B ,Š T )(Š T − S B ) = h(Š T ) − h(S B ) = h(S T ) − h(S B ) (where c(·, ·) is the shock speed introduced in (3.9)), which is precisely equal to the difference of fluxes at the left and the right boundaries.

Conclusion
In this paper, we considered the hyperbolic Buckley-Leverett equation (1.1) and its parabolic counterpart (1.3) for the case when the flux function h is non-monotone. This occurs in gravity driven flows. In the presence of hysteresis, represented by (1.5), the hyperbolic (vanishing capillarity) limit of solutions to (1.3) differs from the classical solution obtained through the equilibrium expression (1.4). In particular, a solution to the Riemann problem (1.1)-(1.2) has a stationary shock at x = 0 if the left (x < 0) and the right (x > 0) states lie in the increasing and decreasing parts of the flux function respectively, or vice versa. The hysteretic states on the left and right become different in this case and thus, different capillary relations are used. Using travelling wave solutions, an admissibility (entropy) condition (3.24) is derived for stationary shocks. The condition states that the flux function h and the pressure corresponding to the hysteretic state of the system, remain continuous across the shock. It is then used to derive all admissible shocks. They are listed in Table 1. The shocks are classified into two categories. The first (Case A) connects imbibition states to drainage states, and the second (Case B) has one of the states undetermined.
Depending on the values of the Riemann data with respect to characteristic points that are easily computable a-priori from h and the capillary curves, two possibilities are identified (corresponding to Case A and Case B) when the stationary shocks occur. The vanishing capillarity solutions for these cases are given by (4.8) and (4.10). Interestingly, the solution (4.10) remains frozen in time in one of the halves and thus, differs significantly from the classical solution. To our knowledge, this is a novel observation. The predictions were validated using numerical experiments.