Fronts in two‐phase porous media flow problems: The effects of hysteresis and dynamic capillarity

In this work, we study the behavior of saturation fronts for two‐phase flow through a long homogeneous porous column . In particular, the model includes hysteresis and dynamic effects in the capillary pressure and hysteresis in the permeabilities. The analysis uses traveling wave approximation. Entropy solutions are derived for Riemann problems that are arising in this context. These solutions belong to a much broader class compared to the standard Oleinik solutions, where hysteresis and dynamic effects are neglected. The relevant cases are examined and the corresponding solutions are categorized. They include nonmonotone profiles, multiple shocks, and self‐developing stable saturation plateaus. Numerical results are presented that illustrate the mathematical analysis. Finally, we discuss the implication of our findings in the context of available experimental results.


INTRODUCTION
Modeling of two-phase flow through the subsurface is important for many practical applications, from groundwater modeling and oil and gas recovery to CO 2 sequestration. For this purpose, the mass balance equations are used which read in the absence of source terms 1,2 as follows: ( ) where = denotes the nonwetting phase and = the wetting phase. Furthermore, is the porosity, and the saturation and density of the phases. The phase-velocities are given by the Darcy's law, 1 Here [m 2 ] is the absolute permeability of the porous medium, [Pa ⋅ s] the viscosity, and the relative permeability of each phase. Moreover, [Pa], [m∕s 2 ], and̂stand for the phase pressure, the gravitational acceleration, and the unit vector along gravity, respectively. Observe that the system (1)- (2) is not closed as there are more unknowns than equations, that is, , , and . Hence, one needs to take certain assumptions. Assuming incompressibility results in being constant. Moreover, by definition Commonly it is assumed that the relative permeabilities, as well as the phase pressure difference, are functions of the saturation of the wetting phase, 1 The function ∶ (0, 1] → ℝ + is referred to as the capillary pressure function. System (1)-(4) reduces to the hyperbolic Buckley-Leverett equation if this term is neglected, that is, ≡ 0. The model given by (1)-(4) works well under close to equilibrium conditions and when flow reversal does not take place. However, some more general cases cannot be explained by this model. One of the first evidences of deviation from the standard model was reported in the 1931 paper by Richards 3 where he concluded that the capillary pressure term is hysteretic in nature. Capillary hysteresis refers to the phenomenon that measured for a wetting phase infiltration process follows a curve, denoted here by ( ) ( ), which differs from measured for a drainage process, denoted by ( ) ( ). If the process changes from infiltration to drainage or vice versa, then the follows scanning curves that are intermediate to ( ) ( ) and ( ) ( ). 4 This is shown in detail in Figure 1 (left). Since then, hysteresis has been studied experimentally, [5][6][7] analytically, [8][9][10][11][12] and numerically. [12][13][14][15][16] Variety of models have been proposed to incorporate the effects of hysteresis, such as independent and dependent domain models [17][18][19] and interfacial area models. [20][21][22][23] A comprehensive study of these models can be found in Ref. 24. Using thermodynamically constrained averaging theory (TCAT), 25,26 one can eliminate hysteresis in capillary pressure altogether by introducing interfacial area and Euler characteristics as additional unknowns. A mathematical study of such models is undoubtedly interesting. However, they require additional constitutive equations. Thus, their analysis is beyond the scope of the present work. In this paper, we will use the play-type hysteresis model 4,24 that approximates scanning curves as constant saturation lines. Such models are generally implemented in classical porous media simulators. Well-posedness results for play-type hysteresis model are found in Refs. 8,9,11,12,27. It has a physical basis 4,28 and it can be extended to depict the realistic cases accurately. 24 A similar hysteretic behavior is observed for the relative permeabilities too, although to a lesser extent. Hysteresis of the nonwetting phase relative permeability in the two-phase case (oil and water for example) is reported in Refs. [29][30][31]. The wetting permeability also exhibits hysteresis 6,32 but the effect is less pronounced (see Figure 1 (right)).
Another effect that cannot be explained by the standard model is the occurrences of overshoots. More precisely, in infiltration experiments through initially low saturated soils, it is observed that if the flow rate is large enough then the saturation at an interior point is larger than that on the boundary even in the absence of internal sources. [33][34][35] This cannot be explained by a second-order model such as (1)- (4). [36][37][38] Hence, based on thermodynamic considerations the dynamic capillary model was proposed in Ref. 39. Since then the dynamic capillary term has been measured experimentally 40,41 and it was used successfully to explain overshoots. 9,11,12,[42][43][44][45] Also the well-posedness of the dynamic capillarity model has been proved [46][47][48][49] and numerical methods have been investigated. [50][51][52][53][54] In this paper, we are interested in studying how the flow behavior is influenced if one considers the nonequilibrium effects, that is, hysteresis and dynamic capillarity. For this purpose, we study the system in a one-dimensional setting. The one-dimensional case is relevant when one spatial direction is dominant; it approximates flow through viscous fingers 12,55,56 and it can explain results from the standard experimental setting shown in Figure 2. [33][34][35] In this study, the behavior of the fronts is investigated by traveling wave (TW) solutions. The TW solutions can approximate the saturation and pressure profiles in infiltration experiments through a long column, and the existence conditions of the TWs act as the entropy conditions for the corresponding hyperbolic model when the viscous terms are disregarded. For the unsaturated case ( = 0), TW solutions with dynamic effect were analyzed in Ref. 42 For the two-phase case, it was shown rigorously in Refs. 43, 44, 57 that nonmonotone TWs and nonstandard entropy solutions are existing if one includes dynamic capillarity effect. Similar analysis but for higher-order viscous terms containing spatial derivatives were performed in Refs. 58,59. The existence of TW solutions for the unsaturated case when dynamic capillarity and capillary hysteresis are present was proved in Refs. 9, 45 and criteria for nonmonotonicity and reaching full saturation were stated. It is evidenced in Refs. 15, 60 that hysteresis can explain stable saturation plateaus but it cannot initiate them. In Refs. 10, 12, it is shown how both hysteresis and dynamic capillarity are F I G U R E 2 Setup of an infiltration experiment. Water is injected at a constant rate at the inlet of a column having height . The main axis of the column is orientated along the gravity vector required to explain the growth of viscous fingers. The entropy conditions for Buckley-Leverett equation considering hysteresis in only permeability were derived in Refs. 61-63. However, hysteresis and nonlinearities were not included in the viscous term. This is taken into consideration in Ref. 64 where the authors add a dynamic term to model permeability hysteresis, while disregarding hysteresis and dynamic effects in capillary pressure. The behavior of TW for a nonmonotone flux function in the presence of a third-order term was described in Ref. 65.
In our current work, we build upon 9,42-45 to describe the behavior of fronts in the two-phase case when dynamic capillarity and both type of hysteresis are included. The models that are used in our analysis are introduced in Section 2. Section 3 discusses the existence of TWs when hysteresis and dynamic effects are included in the capillary pressure but not in the permeabilities. Entropy conditions are derived and they reveal that there can be nonclassical shocks. In Section 4, the analysis is extended to include hysteresis in the permeabilities. This makes self-developing stable saturation plateaus and a broader class of entropy solutions possible. Section 5 presents numerical results that support our analytical findings. Finally, we make some concluding remarks in Section 6 and compare the results with experiments.

MATHEMATICAL MODEL
This section is dedicated to the formulation of a mathematical model that can be used to describe an infiltration process of a fluid into a homogeneous porous column. An example for such an infiltration process is the injection of water into a dry sand column (see Figure 2).

Governing equations
Here we consider the one-dimensional situation where the flow problem is defined on an interval (0, ). This simplification is justified by the fact that the walls of the porous column, in which the fluids are injected are impermeable and that saturation is in general almost constant across the section area of the column. The axis is pointing in the direction of gravity. The medium is assumed to be homogeneous and the fluids are incompressible. Under these constraints, (1)-(2) simplify to where and are denoting the time and space variables, respectively. To further simplify the model, after summing (5) for the two phases and using (3) we observe that the total velocity is constant in space. In addition to that, we assume that is also constant in time, which occurs, for example, if a constant influx (injection rate) is prescribed at the inlet = 0. This gives = = 0, or ( , ) ≡ for ∈ (0, ) and > 0.
From (5)- (7), one finds At this stage, we define the fractional flow function and the function Substituting these relations and definitions into (5) for = yields the transport equation for the wetting phase where we used the notation ∶= and ∶= − .
Note that and ℎ are functions of and possibly of , as shown below.

Modeling hysteresis and dynamic capillarity
To incorporate hysteresis and dynamic capillarity in the model, one needs to extend capillary pressure and relative permeability given in the closure relationship (4).

Capillary pressure
The following expression is used to extend the capillary pressure: where sign(⋅) denotes the multivalued signum graph (see Refs. 4,11,39). The second and third terms in the right-hand side of (13) describe, respectively, capillary hysteresis 4 and dynamic capillarity. 39 Furthermore, ≥ 0 denotes the dynamic capillary coefficient. It models relaxation or damping in the capillary pressure. Although in practice may depend on , 34 Here, and later in this paper, a prime denotes differentiation with respect to the argument. In the absence of dynamic effects, that is, = 0, expression (13) implies This is precisely what is seen from water infiltration/drainage experiments. 5 When = 0, is between ( ) ( ) and ( ) ( ). For this reason, the hysteresis described by (13) is called play-type hysteresis: that is, the scanning curves between ( ) ( ) and ( ) ( ) are vertical.
Before discussing the case > 0, we introduce for convenience the sets and the strip  =  ( ) ∪  ∪  ( ) = {0 < ≤ 1}. In Ref. 45, it is shown that pressure expression (13) can be written as

Relative permeability
To make the effect of hysteresis explicit in the relative permeabilities, we need to incorporate a dependence on both and . This dependence should satisfy Here, ( ) , ( ) ∶ [0, 1] → ℝ are the infiltration and drainage relative permeabilities obtained from experiments. 6,[29][30][31][32] In line with the experimental outcomes, we assume here for ∈ { , }, Note the reverse ordering in and when switching from infiltration to drainage. This is demonstrated experimentally in Refs. 29, 30, 32, see also Figure 1.
In Ref. 66, a play-type approach has been proposed to model where However, this model is ill-posed in the unregularized case as for = 0 the relative permeabilities are undetermined, that is, the relative permeabilities have no equation to determine them when = 0. This is different for the capillary pressure (13) because satisfies Equation (11) as well. With the permeabilities, we take an approach inspired by Refs. 61-63. Here, inherited from the capillary pressure, the hysteresis is of the play-type as well, but now depending on and , rather than on and . We propose the following model: for Here,̄∶  → [0, ∞) is a given function that satisfies Observe that this implies̄( , ( ) ( )) = ( ) ( ) for ∈ { , }. For the moment we leave the choice of̄unspecified except for properties (P4), as it neither influences the entropy conditions nor the critical values introduced afterward.
Using (19) and (9),(10), the nonlinearities and ℎ are expressed in terms of and as well: From (P2)-(P4) we deduce for and ℎ: Observe that in general no ordering holds between ℎ ( ) and ℎ ( ) . Typical curves for ( ) and ℎ ( ) are shown in Figure 3. The Equations (11), (13), and (21) are the complete set of equations for our model.

Dimensionless formulation
The closure relation (13) becomes Now choosing = = 2 2 , the Leverett scaling for gives This choice leaves us with a characteristic dynamic coefficient that is independent of the length scale of the problem. This is precisely the scaling used in Refs. 43, 44, 68 that is consistent with the hyperbolic limit. Realistic values of dimensional and scaled quantities are given in Ref. 69. Dropping thẽsign from the notation, we are left with the dimensionless system This system can be seen as a regularization of the hyperbolic Buckley-Leverett equation with gravity.
Here, the regularization involves hysteresis and dynamic capillarity. Compared to the usual secondorder parabolic regularization, yielding shocks that satisfy the standard Oleinik conditions, 70 different (nonparabolic) regularizations may yield shocks that violate these conditions, see, for example, Refs. 43, 71. Such shocks are called nonclassical. One of the main issues that concerns this paper is to show the existence of nonclassical shocks originating from System (). To this end, we proceed as in Ref. 43 and study the existence of TW solutions of () that connect a left state to a right state in the presence of both hysteresis and dynamic capillarity. TWs for the model with only dynamic capillarity are analyzed in Refs. 44, 57. For the case of unsaturated flow, that is, Richards equation with a convex flux function, existence and qualitative properties of TWs are considered in detail in Refs. 9,42,45. For the purpose of TWs, we consider System () in the domain −∞ < < ∞. Then the capillary number can be removed from the problem by the scaling ∶= ∕ and ∶= ∕ .
This yields the independent formulation with −∞ < < ∞ and > 0. This is the starting point for the TW analysis. These properties of and will be used when discussing the different cases of TWs.

TW formulation
Having derived the nondimensional hysteretic two-phase flow System (), we investigate under which conditions TW solutions exist. These are solutions of the form where and are the wave profiles of saturation and pressure and ∈ ℝ the wave-speed. We seek TWs that satisfy where corresponds to an "initial" saturation and to the injected saturation. The choice of ′ (±∞) = 0 ensures that the diffusive flux vanishes at = ±∞. Substituting (28) into (26), and integrating (26) one obtains where ∈ ℝ and is a constant of integration. As was shown in Ref. 9 for the Richards equation, (28) and (29) do not automatically guarantee the existence of lim →±∞ ( ). However, if (±∞) is well-defined then (29b) and the existence of (±∞) We show later that , interpreted as the initial pressure, can sometimes be chosen independently, whereas, , when existing, is always fixed by the choice of , , and . Following the steps in Refs. 9, 44, 45, we obtain the Rankine-Hugoniot condition for wave-speed , that is, With this, system (29) can be rewritten as a dynamical system, where Note that when is nonmonotone (eg, > in Figure 2), the wave-speed can be positive or negative depending on the values of and . We study all possible solutions of system (TW) for > 0. They serve as viscous profiles of admissible shocks of the limiting Buckley-Leverett equation. Existence conditions for solutions of (TW) act as admissibility/entropy conditions for the corresponding shocks.
The solutions of (TW) are investigated under two different scenarios.
A: No hysteresis in relative permeabilities, that is, ( ) = ( ) for ∈ { , ℎ}. Furthermore, is sufficiently small so that satisfies properties stated for ( ) in (P5). For as in Remark 2 this is satisfied if ≤ . A third scenario where is large so that is nonmonotone is discussed briefly at the end of Section 3.

NO RELATIVE PERMEABILIT Y HYSTERESIS AND SMALL (SCENARIO A)
In the absence of relative permeability hysteresis, the functions , ℎ, , and  depend on only. We explicitly state the properties of as a result of (P5), (P6), and Remark 2.

Preliminaries
For the purpose of readability, we moved the proofs of most details in Appendix A. Throughout this section, we restrict ourselves to relatively small values of . Specifically, we assume First let us take ≤ , where is the saturation at which ′ ( ) = To each ∈ [̃,̄] corresponds a unique ∈ [̄, 1] such that ( , ( )) is the third intersection point between the graph of and the chord through ( , ( )) and ( , ( )) (see Figure 4 (right)). This Values of and are also shown. Note that ( 1 ) and ( 4 ) do not exist Later in this section, a second function = ( ) is introduced as one of the roots of the equation Here, ( ; , ) is the expanded notation of  from (32) for the independent case: A typical sketch of ( ; , ) for different values of is shown in Figure 5. Note that Since, when ≠̃, as ↗ 1, we have for most practical applications The functions , , and̄(assuming (39)) and the definitions of * and * for < This is the case for Brooks-Corey permeabilities with ≥ 2 (see Remark 2).
Remark 3. We also note that in the context of this section the wave-speed (30) reduces to From assumption A1, it follows that there is a one-to-one correspondence between and ∈ [ ,̄]. Writing = ( ), we have the following properties for ( ): Returning to Equation (35), we note that = is the trivial solution. Properties (37) and (38) imply the existence of a second nontrivial and increasing solution = ( ) for ≥ . It satisfies the following properties. Proposition 1. Assume either (39) or ∈ ( , ). Let be the increasing (unique) solution of (35). Figure 6).
Remark 4. For simplicity, we assume (39) for the remaining discussion. This guarantees the existence of a ( * , * ) pair. The methods presented in this paper can also be applied to analyze the case when ( ) and ( ) are not intersecting. The results are briefly discussed in Section 3.2.
Next we turn to system (TW) where, for the time being, we take The − phase plane and the direction of orbits for Scenario A with̃≤ ≤̄. The regions ,  ( ) ,  ( ) and the equilibrium lines are marked Since (TW) is autonomous, it is convenient to represent solutions as orbits in the ( , )-plane, or rather, in the strip {( , ) ∶ 0 ≤ ≤ 1, ∈ ℝ}. Moreover, orbits are same for any shift in the independent variable . Therefore, we may set without loss of generality (see Refs. 9, 45), Equilibrium points of (TW) that are of particular interest are As in Refs. 9, 45, we deduce from (TW) that should satisfy Using techniques from, 9,45 one can show that initial value problem (43), (42a) has a unique local solution ( ; , ) that satisfies (42b).
Remark 5. Conversely one recovers the orbit ( ( ), ( )) by substituting into (31). Using (41)  Rewriting (43) as we find recalling (P1) that Integrating this inequality from to gives the lower bound With satisfying (40), properties (37)-(39) and Proposition 1 imply and Observe that, depending on , , and , the interval where ( ) < ( ) ( ) is either ( , 1] if ( ) and ( ) ( ) do not intersect, or ( , ) with ≤ 1 in case there is an intersection at = . In the latter case, it follows immediately from (44) that we must have, Hence, if the orbit exits through the capillary pressure curve ( ) , it can only do so at points where  ≥ 0. From the discussion above, one defines which is the upper limit of the interval on which exists. Then we have In Figure 8, we sketch the behavior of ( ; , ) in  ( ) . The existence of orbits as in Figure 8 (left) is a direct consequence of the behavior of the lower bound (45). Orbits as in Figure 8 (right) need more attention since the case ( , ) < ( ), represented by 3 , remains to be discussed. We make the behavior as sketched in Figure 8 (right) precise in a number of steps.
Next, we give a general monotonicity result.
So far we have shown monotonicity of the orbits in  ( ) . The next results give the continuous dependence on the parameter and . This is addressed in the following results.
Proposition 5 (Continuous dependence of ). Let = ( ) − . In the context of Figure 4 and with Φ defined in (45), we have Corollary 1 (Continuous dependence of ). Let 0 > 0 and 0 be fixed such that Ordering of the orbits in the -phase plane for * < <̄and typical values of : The behavior of the orbits in the -plane for the same values Now we are in a position to describe how behaves for different combinations of and .
The proof of this result is given in proposition 2.1 of Ref. 45. For later purposes, we introduce here some definitions that are related to Proposition 6. Then, which directly gives ( , ) = . Thus, > 0 is defined as Finally, using [Ref. 9, proposition 4.2(b)] one concludes that ( ) < ∞.

Problem (TW) with
< ≤W e investigate the existence of an orbit connecting ( , ( ) ( )) and the segment . Defining . We omit the details here. In Theorem 1, we have taken < without loss of generality. In the > case, the roles of the equilibrium points and are reversed. The typical behavior of the and the profiles with respect to is given in Figure 11. Both and are monotone for < , whereas for < < they have finite number of local extrema and (+∞) = ( ) ( ). For > , F I G U R E 1 1 Typical behavior of ( ) (left) and ( ) (right) for different values of . Here profiles for three different values are plotted satisfying 0 < 1 < < 2 < < 3 . The = 0 coordinate is fixed by (41) F I G U R E 1 2 The sets , , and  and the functionš( ),̂( ) has infinitely many decaying local extrema, whereas has no limit. In particular, each maximum corresponds to a saturation overshoot. On the other hand, the oscillations in become wider, in line with the assumption lim →∞ ′ ( ) = 0. In this case, the segment becomes an -limit set of the orbit.
Turning to the case, ∈ ( * ,̄), we define the two functions which will be used extensively below. Observe thať( ) is a strictly decreasing function whereaŝ( ) is a strictly increasing function for >̄and̂( ) =̌( ) =̄for ≤̄. This is sketched in Figure 12. Numerically computeď ( ) and̂( ) functions are shown in Figure 20. With this in mind, we define the following sets: Observe that if <̄then only regions  and  are relevant. With defined in (A1), for ( , ) ∈  one has We discuss the remaining situations, ( , ) ∈  and ( , ) ∈  in the next section.

Since
>̌( ), a TW cannot connect and . However, a different class of waves is possible when ( , ) ∈ .

Proposition 9. For a fixed
∈ (0, ) and ( , ) ∈ , consider the system For this system, an orbit ( , ) exists that connectŝ( Proof. Upon inspection of the eigendirections for the system (52) around the equilibrium point̂( ) one concludes that there is indeed an orbit ( , ) that connects tô( ) as → −∞ from the set  ( ) defined in (15). Moreover, from the direction of the orbits in this case, as shown in Figure 13 (left), it is apparent that after leavinĝ( ) , decreases monotonically till the orbit either hits the curve = ( ) ( ) for some ≤ or exits { > } through the line = . We prove that it is not possible for the orbit to escape through = . To show this, consider the orbit ( 1 , 1 ) that satisfies the original (TW) equations and enters  ( ) from̂( ) . We show that this orbit cannot cross the line = .
The divergence argument presented in Refs. 9, 42, 45 is used for this purpose. To elaborate, assume that ( 1 , 1 ) intersects the line = at . Consider the region Ω, enclosed by the segments , , the orbit ( 1 , 1 ) and the orbit ( 2 , 2 ) that satisfies (TW) and connects tô( ) (see Figure 13 (right)). Introducing the vector-valued function ⃖⃖ ⃗ ( , ) = ( 1  ( , ), ( ; , )) and deduces from (16), This gives a contradiction when the divergence theorem is applied to ⃖⃖ ⃗ in the domain Ω: the integral of ⃖⃖ ⃗ over Ω is nonnegative whereas ∫ Ω div ⃖⃖ ⃗ < 0 from (P1) and Figure 13 (right). Hence, the orbit The wave-speed corresponding to the orbit ( , ) satisfies being the speed of both ( 1 , 1 ) and ( 2 , 2 ) waves. Hence, by the continuity of the orbits with respect to , as shown in Proposition 4, it is evident that ( , ) intersects ( ) ( ) for some > . Observe that, if ( , ) ∈  then TWs do not exist between and̂( ) since both are in the concave part of with >̂( ). Thus we have exhausted all the possibilities of connecting and with Theorem 1 and Proposition 9.

Entropy solutions to hyperbolic conservation laws
Under the conditions of Scenario A, we consider the Riemann problem In the context of the viscous model discussed in this paper, we consider the Buckley-Leverett Equation (54a) as the limit of System (25) for ↘ 0. As a consequence, we only take into account those shock solutions of (54a) that have a viscous profile in the form of a traveling wave satisfying (TW). Such shocks are called admissible because they arise as the → 0 limit of TWs. In this sense, the entropy condition for shocks satisfying (54a) are equivalent to existence conditions for traveling waves satisfying (TW). This may lead to nonclassical shocks violating the well-known Oleinik entropy conditions, see, for example, Ref. 43.
Here, we assume which is more general compared to (33) where the additional constraint of < was imposed. This generalization is possible since > simply implies that the sets ,  are empty. Our analysis can also be applied to derive the entropy conditions for the case > , however, for simplicity we restrict our discussion to (55).

( , ) ∈ 
As in the usual Buckley-Leverett case (ie, without dynamic capillarity and hysteresis in the regularized models) the solution is given by Here, the shock satisfies the classical Oleinik condition.

( , ) ∈ 
In this case, the admissible solution is composed of two shocks: an infiltration shock from tô( ), followed by a drainage shock from̂( ) to .
Note that this solution is different from the Oleinik entropy solution. 70 Moreover, the shock connecting and̂( ) is undercompressive 71 meaning that it violates the Lax entropy condition. This is because > ′ (̂( )) for this shock. The trailing shock connectinĝ( ) to is however classical.

( , ) ∈ 
The solution in this case violates again the Oleinik entropy condition. It consists of an infiltration shock from tô( ) followed by a rarefaction wave from̂( ) to , with (⋅) satisfying Since is concave for ∈ [̂( ), ], ′ is monotone implying that (⋅) is well-defined. We observe that in the last two cases the solution features a plateau-like region. This plateau appears and grows in time since the speeds of the drainage shock and of the endpoint of the rarefaction wave are lesser than the speed of the infiltration shock. Interestingly, the saturation of the plateau only depends on ( ) and not on ( ) . To be more specific, although the viscous profile consisting of a TW wave connectingâ nd depends on ( ) , the shock solution resulting from it, in the hyperbolic limit, does not. However, the role of the drainage curve in the entropy solutions become evident in Scenario B, which is discussed in the next section.
In the absence of hysteresis and for linear higher-order terms, which correspond to constant and linear -dependence, in [Ref. 43, section 6] it is proved that the nonstandard entropy conditions discussed here are entropy dissipative for the entropy ( ) = 1 2 2 . However, such an analysis is beyond the scope of this paper. The solution profiles for the Riemann problem are shown in Figure 14.

Extension to the nonmonotone case
The analysis so far can be extended to the case where is large resulting in being nonmonotone. If ∈ (0, 1) is the saturation where ( ) attains its maximum (see Remark 2 and Figure 3), then the results obtained so far cover the case when and * are below . However, if > then the TW study has to be conducted also from an perspective, not only from the one. In this scenario, since fronts having negative speeds and thus moving toward become possible, one has to consider the functionŝ( ),̌( ) for a fixed , similar tô( ),̌( ) from Definition 1 for fixed . Due to the symmetry in the behavior of the fronts approaching , respectively, , some of the results obtained so far extend straightforwardly to the nonmonotone case. However, a detailed analysis is much more involved and therefore left for future research because of the following two reasons: 1. Depending on the relative positions of ,̂, , and̂, there are many subcases to consider. In this case, up to three shocks are possible, traveling both forward and backward. Which of these shocks are admissible and how they are connected requires further analysis.
2. For a nonmonotone , when considering the hyperbolic limit in the absence of hysteresis or dynamic effects, the entropy solutions may include rarefaction waves with endpoints moving in opposite directions, forward and backward. When capillary hysteresis is included, preliminary numerical results have provided solutions incorporating two rarefaction waves, one with endpoints traveling backward and another one with endpoints traveling forward, and a stationary shock at = 0. Such solutions still need to be analyzed further.
In this scenario, can be taken in the entire interval (0,1) and can be chosen independently as long as ( , ) ∈ , that is, This is different from Scenario A where is restricted to the interval (0, ) and is fixed to = ( ) ( ). We first introduce some notation. Regarding the choice of , we have the following.

Proposition 10. Let and be as in Case (i) or Case (ii). Then any solution of (TW) that connects and can only exist if
Proof. Since is an equilibrium point,  ( , ) = 0, which implies that ∈ [ ( ) ( ), ( ) ( )]. The directions of the orbits for in this interval are displayed in Figure 15 (right). We proceed by introducing the set It corresponds to the black dotted curve in Figure 15 (60). Concerning ( ) , has either zero, one or two intersection points (see Figure 16 (left)). In the latter case, as before,  0 contains one or two vertical half-lines as shown in the (right) plot of Figure 16. However, this aspect plays no major role in the analysis below.
Every point in the set  0 ∩  is an equilibrium point. However, all points in the set  0 ∩ int() (the interior of  being referred to as int() here) are unstable and as follows from Figure 15 (right), no orbit can reach these points as → ∞. This eliminates all other possibilities to reach as → ∞ except for = ( ) ( ) and = ( ) ( ). ■ We now consider the two cases separately.

Case (ii):̄≤ < and stability of plateaus
The counterpart of Proposition 11 for Case (ii) is (see also Figure 17 Finally, we investigate a special case related to the development of stable saturation plateaus in infiltration experiments. For ∈ (0, 1), and ∈ ( , 1) a stable plateau is formed when an infiltration wave, from to ∈ ( , 1), followed by a drainage wave, from to , both have the same speed resulting in the width of the plateau to remain constant. This is different from the plateaus described in (57) where the speeds of the infiltration and the drainage fronts are necessarily different. The existence of stable saturation plateaus has been widely studied experimentally 29,33,35 and numerically. 15,60 Although results regarding stability of the plateau are available, 15,60 the mechanism behind its development is still not well understood. Here, we give an example where our analysis predicts that such a plateau will develop. Specifically, it occurs when > * ( ) and a direct monotone orbit from to ( , ( ) ( )) is no longer possible. This is verified numerically in Section 5.2. The proof follows directly from Propositions 11 and 12.

Entropy solutions
We can now discuss the entropy solutions of the Riemann problem (54) under the assumptions of Scenario B. To be more specific, we give a selection criteria for the solutions of the system with We view (63) as the limit of () when the capillary effects vanish. However, hysteresis is still present in the model. Note that still plays a role in determining the entropy solution despite being absent in (63). This is similar to what we saw in Section 3. However, the focus here being hysteresis in permeability and capillary pressure, for a fixed ∈ (0, 1) we take Observe that (65) does not provide a void interval for . To see this, note that * ( ) is defined similar to in Proposition 6 and thus, it satisfies the inequality in (49), that is, it has the positive quantitȳ as its lower bound. Although̄in Proposition 6 actually depends on , one sees from (48) that the values of̄are bounded away from 0 uniformly with respect to . Hence, * ( ) is also bounded uniformly away from 0. Similar argument holds for * ( ).
We now consider the cases > and < separately.

NUMERICAL RESULTS
For the numerical experiments, we solve () (System (26)) in a domain ( , ), where < 0 and > 0. As an initial condition for the saturation variable, we choose a smooth and monotone approximation of the Riemann data: Here, is a smoothing parameter, denotes the saturation induced by a certain injection rate and is the initial saturation within the porous medium. To model the capillary pressure, a van Genuchten parameterization is considered, that is, In the remainder of this section, we use the following parameter set: Λ = 3.5, = 0.92, Λ = 7, and = 0.9. To solve () numerically, for ∈ ℕ ∪ {0} and 0 = 0, we solve within the time step with respect to the pressure variable . For a given , this is a nonlinear elliptic problem and to solve it, a linear iterative scheme is employed which is referred to as the L-scheme in literature [73][74][75] : Here, denotes the pressure at the th iteration and 0 = ( , ). On closer examination, the L-scheme corresponds to a linearization of the nonlinear problem, since for each iteration a linear equation in the unknown pressure variable is solved. For Scenario A the parameter is set to = 1 to ensure convergence of the L-scheme 45,73 and for Scenario B the modified variant of the L-scheme is used 9,75 to speed up the convergence, since in this scenario the stiffness matrix has to be recomputed in every iteration. A standard cell centered finite volume scheme is considered for discretizing the linearized elliptic problem in space. Having the pressure variable and the saturation variable for = at hand, we update the saturation as follows: .

Numerical results for Scenario A
First we illustrate the theoretical findings of Scenario A. The boundary conditions with respect to the pressure variable are of Neumann type at = and of Dirichlet type at = : The boundaries of the domain are given by: = −10 and = 500. Because we do not include hysteresis in the relative permeabilites, the flux function depends only on and is determined by: The numerical results presented in this subsection are related to = end = 300. For the parameters of the initial condition, we take: = 0.1, = 0.4, and = 1.
Moreover, the curves for and are determined (see Figure 18). Observe that, from our choice, < and ∈ ( , * ]. Next, the characteristic -values for drainage and imbibition are computed. Using (50) and given parameters, we obtain: and study the resulting -orbits. Considering Figure 19, it can be observed that for < monotone The parameter choice considered so far, corresponds to the solution class  (see (51)), whose entropy solution consists of a single shock without any saturation overshoots (see Figure 21 (top)). However, there are two further solution classes,  and  (see (51)), arising in the context of Scenario A, represented by entropy solutions (57) and (58). In case of solution class , the entropy solution is given by saturation plateau that is formed by an infiltration wave followed by a drainage wave. The saturation at plateau level is denoted bŷ( ). For solution class , the entropy solution exhibits a rarefaction wave connecting witĥ( ), which is connected to by a shock. To observe these cases numerically, we compute thê( ) anď( ) curves introduced in Definition 1 (see Figure 20). In the figure, we fix = 1 and vary so that the pairs ( , ) belong to one of the sets , , and . The results are shown in Figure 21 with the (left) plot showing the variation of with , and the (right) plot showing the profiles in the -phase plane. The curves corresponding to Set  show a direct TW wave connecting and = 0.35. Some oscillatory behavior around can be observed since is comparatively large, however, the existence of a single TW between and implies that these states are connectable by an admissible shock in the hyperbolic limit. Next, choosing = 0.55, ( , ) lies in Set , and a solution consisting of an infiltration wave followed by a drainage wave is computed in accordance with the theory. The development of the plateau is shown in Figure 22. Again, small oscillations are seen in the drainage wave part which is expected from Corollary 2 since is large. The resulting plateau has saturation 0.7158, whereas, the prediction from Figure 20 iŝ( ) = 0.7254. Finally, for = 0.8, the pair ( , ) belongs to the Set . The numerical solution exhibits a shocklike structure followed by a plateau and they coincide with the infiltration wave of Set  on both plots of Figure 21. Moreover, a rarefaction wave between̂( ) and is detected. Thus, we conclude that the saturation profiles in Figure 21 correspond to the entropy solutions depicted in Figure 14 and the numerical results are in agreement with the theory.

Stability of the TW profiles
Another important issue with regards to our discussion is the stability of the TW solutions. For the dynamic capillarity model, linearized stability of TW profiles was investigated in Refs. 76, 77 for convex flux functions. The TWs were found to be stable under small perturbations. Linear stability of the steady-state solution for the dynamic capillarity model was shown in Ref. 37. However, such a result might not hold when hysteresis is included. For example, in Ref. 38 it was shown that with hysteresis, stable planar fronts are not guaranteed. Moreover, linearized stability results might not extend to large deviations from the TW profile. Nevertheless, in our case, the stability of the TW profiles is indicated by the numerical computations. More precisely, the numerical examples are presenting the solution profiles for the partial differential equations, thus not for the TW. These numerical results are obtained t t t t F I G U R E 2 2 The time evolution of the solution profile for ( , ) ∈  by using a regularized initial data in (70). The solution quickly develops into profiles that are either an approximation of the corresponding TW if the values of and are states that can be connected by a TW for the chosen parameters, or it exhibits a plateau value that approximates very well the left state that can be connected by a TW to . Moreover, this is obtained irrespective of the choice of the regularization parameter in the initial condition (70). Similar behavior was reported in the numerical examples of Refs. 9, 43-45. However, instability of the TWs is observed for large when hysteresis is included in the permeabilities. This requires further attention and we have excluded this possibility in our current analysis and computations.

Numerical results for Scenario B
In case of Scenario B, we choose the following boundary conditions with respect to the pressure variable. As in the previous subsection, they are of Neumann type at = and of Dirichlet type at = : Moreover, the boundaries of the domain are given by: = −10 and = 190. To make matters interesting, contrary to the previous subsection, we do not start with an infiltration state for , but with a drainage state. Due to the fact that we consider hysteresis both in the capillary pressure and relative permeabilities, fractional flow functions are introduced both for infiltration and for drainage. We use , ( ) ( ) =  (67) it is expected that the entropy solution will have a shock from ,1 tō, followed by a rarefaction wave from̄to ,1 . This is exactly what is seen from the viscous profiles obtained numerically (see Figure 23). Similarly, for the second case, since ,2 <̄and is small, we see from Figure 23a viscous solution resembling a drainage shock followed by a rarefaction wave, as predicted in (69). Next, we investigate whether a stable plateau is are located on the same line. This is precisely the condition that the solutions ( ( ) , ( ) ) and ( ( ) , ( ) ) of Proposition 13 satisfy. Drawing a line through the given points for = 0.3 and = 0.5 (see Figure 24), we obtain that a stable plateau should be located at ≈ 0.634. As seen from Figure 24, the orbit in the -plane stabilizes exactly at ≈ 0.634 and all the three points line up. Considering Figure 25, we observe that the saturation plateau is in a transient state in the beginning, but it stabilizes at ≈ 0.634 for longer times, as the speeds of the infiltration and drainage waves match.

FINAL REMARK S AND COMPARISON WITH EXPERIMENTS
In this work, a one-dimensional two-phase flow model has been analyzed for infiltration problems. For simplicity, we have assumed that the medium is homogeneous and a constant total velocity is prescribed at the boundary. Dynamic and hysteretic effects are included in the capillary pressure with transitions between drainage and infiltration processes being modeled by a play-type hysteresis model having vertical scanning curves. Relative permeabilities are modeled as functions of saturation and capillary pressure to make their hysteretic nature explicit. The focus being on TWs, the system of partial differential equations is transformed into a dynamical system. This system is analyzed for two different scenarios, A and B. In Scenario A, the hysteresis appears only in the capillary pressure, and we consider a broad range of dynamic capillarity terms, F I G U R E 2 6 The saturation profiles for water injection experiments taken from Ref. 29. Note the presence of growing saturation plateaus from small to large ones. In Scenario B, hysteresis is included in both the relative permeabilities and in the capillary pressure, whereas the dynamic capillary effects are kept small. For each scenario, the existence of TW solutions is studied. In particular, we show that if the dynamic capillary effects exceed a certain threshold value, the TW profiles become nonmonotonic. Such results complement the analysis in Refs. 9, 42, 45 done for the unsaturated flow case, respectively, in Refs. 43, 44 for two-phase flow but without hysteresis. From practical point of view, the present analysis provides a criterion for the occurrence of overshoots in two-phase infiltration experiments.
Based on the TW analysis, we give admissibility conditions for shock solutions to the hyperbolic limit of the system. Motivated by the hysteretic and dynamic capillarity effects, such solutions do not satisfy the classical entropy condition. This is because the standard entropy solutions to hyperbolic two-phase flow models are obtained as limits of solutions to classical two-phase flow models, thus not including hysteresis and dynamic capillarity. In particular, for the infiltration case of Scenario A, apart from the classical solutions, there can be solutions consisting of (a) an infiltration shock followed by a rarefaction wave having nonmatching speeds, or (b) an infiltration shock followed by a drainage shock resulting in a growing saturation plateau (overshoot) in between. This is similar to the results in Refs. 43, 44. In Scenario B, the entropy solutions are shown to depend also on the initial pressure. In particular, if certain parametric conditions are met, the solutions may include ones featuring a stable saturation plateau between an infiltration front and a drainage front, both traveling with the same velocity. Such solutions are obtained, for example, in Ref. 15, but only after generating the overshoot through a change in the boundary condition. All cases mentioned above have been reproduced by numerical experiments, in which a good resemblance has been observed between the TW results and the long-time behavior of the solutions to the original system of partial differential equations.
From practical point of view, we note that the present analysis can also be used to explain experimental results reported, for example, in Refs. 33,41,55,66. The occurrence of saturation overshoots is predicted theoretically for high enough dynamic capillary effects, namely, of the value in (24). In dimensionless setting, this can be assimilated to an injection rate that is sufficiently large. This is in line with the experimental results in Ref. 33, where the development of plateau-like profiles was observed for high enough injection rates, as shown in fig. 5 of 33 and fig. 5.3 of Ref. 66. Similarly, in the water and oil case, the plateaus are seen to develop and grow in figs. 5-6, 8-9, 18 of Ref. 29 (see Figure 26). This behavior is predicted by the analysis in Section 3. Moreover, fig. 10 of 29 might be presenting the case when the saturation has developed a plateau between two fronts traveling with the same velocity, a situation that is explained by the authors by means of hysteretic effects in the flux functions. Such solutions are investigated numerically in Refs. 15, 60, where it is shown that the plateaus can persist in time but without explaining how they are generated. The results in Section 4 partly support the conclusions there, but also explain the mechanism behind the development of such plateaus. We mention 78 in this regard, where the authors conclude that a similar mechanism must be responsible for observed stable saturation plateaus inside viscous fingers. Proof of Proposition 1. The property that ( ) increases follows directly from (36). Moreover, ( ) > for > , since ∫ ( ; , ) < 0 in this case (see (37)). Observe that, if (39) is satisfied then (̃) < 1. This shows the existence of ( ) in a right neighborhood of =̃. The solution in this case exists up to = * ∈ (̃,̄) where ( ) and ( ) intersect: ( * ) = ( * ) =∶ * .
Since  > 0 in ( , ) for < ≤ ( ), the term in the right of (A4) is negative, yielding a contradiction.  Integrating this equation and using Proposition 4 and Φ from (45) gives the desired inequalities. ■ Proof of Corollary 1. We only demonstrate continuity with respect to . The proof of continuity with respect to follows the same arguments. We therefore take = 0 and drop its dependence from the notation for simplicity. Consider first > 0 and 0 < ( 0 ). Recalling ( ( 0 ); 0 ) = 0, Proposition 5 gives where Φ( ( 0 )) > 0 by (46) and Proposition 3. For any given (small) > 0 and with reference to