Hysteresis and Horizontal Redistribution in Porous Media

It is well known that multiphase flow in porous media exhibits hysteretic behaviour. This is caused by different fluid–fluid behaviour if the flux reverses. For example, for flow of water in unsaturated soils the process of imbibition and drainage behaves differently. In this paper we study a new model for hysteresis that extends the current playtype hysteresis model in which the scanning curves between drainage and imbibition are vertical. In our approach the scanning curves are non-vertical and can be constructed to approximate experimentally observed scanning curves. Furthermore our approach does not require any book-keeping when the flux reverses at some point in space. Specifically, we consider the problem of horizontal redistribution to illustrate the strength of the new model. We show that all cases of redistribution can be handled, including the unconventional flow cases. For an infinite column, our analysis involves a self-similar transformation of the equations. We also present a numerical approach (L-scheme) for the partial differential equations in a finite domain to recover all redistribution cases of the infinite column provided time is not too large.


Introduction
In this paper we consider the flow of two immiscible and incompressible fluids through a homogeneous and isotropic porous medium. It is assumed that the pores of the medium are fully occupied by these fluids. One fluid is the wetting phase and the other one is the nonwetting phase. They are denoted by the subscripts w and n. We disregard the influence of gravity, as we are interested in a horizontal physical system. Further we assume that there are no internal sinks or sources. The corresponding (i.e. macroscopic) equations are well known (Bear 1972;Collins 1976;Helmig 1997) and in dimensionless form they read: (1.1) S n + S w = 1, S n , S w ≥ 0. (1.2) In (1.1), S i , P i and F i denote saturation (assumed to be scaled such that 0 ≤ S i ≤ 1), pressure and volumetric flux of phase i, and k i is the relative permeability of the porous medium with respect to phase i. If the system contains water as the wetting phase and air as the non-wetting phase, then it is generally assumed that P n = P air is constant. Equation (1.1) reduces to Richards equation which is the standard equation that models flow in unsaturated porous domains (Bear 1972;Collins 1976;Helmig 1997): In general, the closure relation between P and S is given in form of an algebraic relationship determined from experiments: where P c : (0, 1] → R is the capillary pressure. There is a vast amount of literature available on the capillary pressure and its interpretation (see references cited above). For properties and closed-form relations we refer to the well-known references (Brooks and Corey 1966;van Genuchten 1980). It is known that multiphase flow in porous media displays hysteretic effects (Richards 1931). The capillary pressure-saturation relationship (henceforth called P c -S w relationship) traces one path while going through an infiltration/imbibition/wetting process and another path while going through a drainage/drying process. This effect can be incorporated into the classical model by replacing P c (S w ) function in (1.4) by two different functions: P im (S w ) for wetting and P dr (S w ) for drying. This works quite well if the porous medium is going through only a wetting or a drying process. But if there is a switch between the two, then the P c -S w curves span the region between the P im and P dr curves in the form of scanning curves (Morrow and Harris 1965). Typical behaviour is shown in Fig. 1. Some of the most common models for porous media hysteresis are listed below.

Similarity Models
The main idea behind this class of models is to express the scanning curves by a closedform relationship that allows the scanning curve to be similar in shape to the imbibition or drainage curves. With the similarity hypothesis, Philip (1964) was able to obtain reasonable approximations for scanning curves. A formally equivalent model was derived by Poulovassilis (1962) based on the independent domain theory. Concepts of domain theory of capillary hysteresis were further used in the models proposed by Mualem (1974Mualem ( , 1984, Mualem and Dagan (1975), Parker et al. (1987) and others [see Viaene et al. (1994), Kool (Zhuang 2017). b Schematic diagram of scanning curves for playtype hysteresis and similarity models (Beliaev and Hassanizadeh 2001;Mualem 1974) (1987)]. The scanning curves obtained in this way are close to the experimental scanning curves. Moreover, this class of models can describe "secondary hysteresis". This refers to the phenomenon that the scanning curve through a point, switching from wetting to drying, is different from the scanning curve through the same point switching from drying to wetting. All other models discussed in this paper do not include this secondary behaviour, because for those models a intermediate point can move back and forth on the same scanning curve. Similarity models can also be used to explain many cases of the horizontal redistribution problem (Heinen and Raats 1999), which is an important benchmark problem for flows in porous media. However, similarity models are not straightforward to apply because in these models the saturation at a point is a function of all previous reversal points (when it switches from infiltration to drainage or vice versa) (Viaene et al. 1994;Kool and Parker 1987). The closedform expressions of the scanning curves actually take in these reversal points as parameters. The order at which the processes (wetting/drying) have gone through plays also a vital role. This leads to book-keeping for each point in space, making the models difficult to handle in practice and in any numerical or analytical approach.

Playtype Hysteresis
In this approach one models scanning curves as vertical lines between P im and P dr (Visintin 2013;Beliaev and Hassanizadeh 2001;Schweizer 2012). To close Eq. (1.1) one replaces (1.4) by an expression of the form where sign(ζ ) is the multivalued signum graph, (1.6) 123 and where P + , P − are defined by: The model based on (1.3) and (1.4) is well posed in the mathematical sense (Lamacz et al. 2011) and can be physically justified by pore-scale (Schweizer 2005) or thermodynamic (Beliaev and Hassanizadeh 2001) arguments. Although this model has the advantage of being simple and local in time (no information on the history of a point is required) the resulting vertical scanning curves do not really resemble the ones from experiments. Moreover, as we show later, the playtype hysteresis model cannot describe all cases of horizontal redistribution.

Interfacial Area Model
Pore-scale simulations have shown (Reeves and Celia 1996) that interfacial areas play an important role in the P c -S w relationship. Motivated by this, a model was proposed in Gray (1990, 1993) and Niessner and Hassanizadeh (2008) based on thermodynamic considerations. The main idea is to introduce the volumetric interfacial area (a wn ) as a new variable, in addition to saturation and pressure, and to assume that a wn is a unique function of saturation and capillary pressure: a wn = a wn (S w , P n − P w ). (1.8) A transport equation for a wn was proposed leading to a new formulation for multiphase flow. The original idea was that this new formulation including interfacial area could account for the full hysteretic nature of the water retention characteristic. However, when analysing the model (Pop et al. 2009) it was shown that for any fixed point x 0 in space, there exists a unique P c -S w curve which satisfies where g is a given smooth function of P c and S w . Thus if S w (x, 0) and P c (x, 0) have only two values, as they do for horizontal redistribution, two P c -S w curves arise. Clearly for general initial conditions, infinitely many P c -S w curves may arise. Hence the concept of primary wetting and drying is not described by the interfacial area model. Hysteresis, in the sense of switching between two curves at a fixed x 0 , could only be included by introducing rate-dependent terms in the coefficients . The purpose of this paper is to introduce a hysteresis model that is based on the playtype approach, having primary drying/wetting curves and in between non-vertical scanning curves. These scanning curves can be chosen in such a way that they are close to experimental data. This model is presented and discussed in Sect. 2. In Sect. 3 the equations and conditions for (horizontal) redistribution are given and in Sect. 4 self-similar solutions are discussed describing all possible redistribution cases using the new model. Then, in Sect. 5 a numerical scheme (L-scheme) for the partial differential equation is presented and computational results are compared with the redistribution cases. Conclusion are given in Sect. 6.

Extended Playtype Hysteresis Model
In this section we will introduce the extended playtype hysteresis model, and discuss its background and some of its properties. Let u denote the capillary pressure, i.e. u := P n − P w . (2.1) We extend the playtype hysteresis by introducing non-vertical scanning curves which become vertical near the saturation end points S w = 0 and S w = 1. To this end we replace ∂ t S w in the sign function in (1.5) by ∂ t (H (S w ) + u), where H : (0, 1) → (0, ∞) is a function chosen in such a way that the corresponding scanning curves have the desired properties. Thus instead of (1.5) we propose where P ± (S w ) and sign(·) are defined as in (1.7) and (1.6), respectively. To better understand what relation (2.2) implies, let us introduce the following sets (see Fig. 2): Definition 1 In the (S w , u) plane we consider the sets In these definitions we considered S w ∈ (0, 1], meaning 0 < S w ≤ 1 to avoid S w = 0, as P im (S w ) and P dr (S w ) become singular at that point. If (S, u) ∈ H, then by (1.7) Hence from (2.2) Definition (1.6) of sign(·) then gives This means that we have both the lower bound conditions are satisfied and so are inequalities (2.5) and (2.6). Moreover, condition (2.7) ensures that scanning curves originating from arbitrary points on ∂H im or ∂H dr remain in H. The question arises how to choose and construct a function H that gives scanning curves close to experimental data and satisfies (2.7) for mathematical consistency. We present a construction that is based on the experiments of Morrow and Harris (1965). Their results for drying and wetting are shown in Fig. 3.
Here the variables are as in Morrow and Harris (1965): the saturation is unscaled (0 < S wr < S w < 1 − S ar ) and the suction, capillary pressure u in this context, is in cm of water. In the construction of H we use the same variables to get a meaningful comparison. We propose for H the form:  Remark 1 For simplicity one could consider dH dS w (S w ) := 1 δ (is constant) for some δ > 0. Then the scanning curves are described by the straight lines: (2.10) By the same argument as above, one would need for each 0 < S w ≤ 1 to ensure inequality (2.7). However this would unrealistically restrict the range of saturations for which the scanning curves remain in H. Beliaev and Hassanizadeh (2001) were the first to derive the extended playtype hysteresis closure relationship. Using thermodynamical arguments they obtained an expression in the inverse form of (2.2): [see Beliaev and Hassanizadeh (2001), eq. 35]. In our notation it translates into: where C > 0 is a constant and P c 0 is a reference curve between P im and P dr that intersects all scanning curves. For example P c 0 could be P + . The authors argued, using experimental data,

123
that the term C P c 0 (S w ) has roughly 10% contribution in expression (2.12). Also they pointed out that C → 0 corresponds to the playtype hysteresis model. As to put (2.12) in the form of (2.2). Note however that the reference pressure P c 0 does not necessarily satisfy condition (2.7).
Extension of playtype hysteresis model by inclusion of non-vertical scanning curves has been hypothesized also in the context of numerical analysis. One of the earliest example of this is found in Hanks et al. (1969) where a numerical method has been proposed for one dimensional hysteretic columns. In Brokate et al. (2012), authors use non-vertical approximations to vertical scanning curves as the original playtype hysteresis model poses difficulties for convergence. Furthermore, in Lamacz et al. (2011) and in van Duijn et al. (2018), the sign function has been regularized for mathematical analysis, making it resemble non-vertical scanning curves.
Remark 2 Because of the appearance of the term ∂ t u in (2.2), the extended playtype hysteresis model requires initial conditions both in capillary pressure and saturation. This is in contrast to playtype hysteresis where only an initial condition in saturation is required.

Redistribution
Following Philip (1991) , Raats and van Duijn (1995) and Pop et al. (2009) we make the effect of hysteresis explicit in cases of horizontal redistribution. We show that the extended playtype hysteresis as described by (2.2), with H (S w ) satisfying (2.7), covers all these cases. This is one of the main purposes of the paper.

General Set Up
Consider a horizontal porous column of infinite extent, directed along the x-axis. In the column the flow is one dimensional and describes the redistribution of fluids. We shall restrict ourselves to the case of flow in an unsaturated porous media, with water as the wetting phase and air as the non-wetting phase. Setting P n = P air = 0, the variables to be determined are the water saturation S w and the suction u = −P w (same as capillary pressure in this context). For brevity we drop the subscript w from the notation. Let R − = {x < 0} and R + = {x > 0} denote, respectively, the left and right half of the column. The governing equations are Richards equation (1.3) along with closure relation (2.2): The halves R − and R + have constant, but different, initial saturation and suction at t = 0. We impose Throughout this paper we assume that the functions k, P im , P dr and H are smooth and satisfy the structural properties: For the initial conditions (3.4) and (3.5) to be consistent with expression (3.3) we must impose, (3.6) In this and later sections we use the notation, Thus in terms of Definitions 1 and 2, condition (3.6) reads: Mass conservation and momentum conservation requires the flux F and the suction u to be continuous. However, the saturation maybe discontinuous across x = 0. Thus the strategy is to solve (3.1)-(3.3) separately for R − and for R + , subject to (3.4) and (3.5) and then to match possible solutions so that flux and suction are continuous at x = 0. Such solutions will be either in the wetting state or in the drying state in R − and R + . Thus in addition to (3.1)-(3.5) we shall explicitly use for each t > 0. Here we use the notation f (0 ± ) = lim x↑↓0 f (x).

Possible Initial Conditions
Let E r ∈ ∂H dr , E l ∈ ∂H im and let u l > u r . Since this implies S l < S r one expects that water flows from the wet half column to the dry half column. This is called "conventional flow". It is described by Philip in his classical paper Philip (1991). He found that in this case the right half column (R + ) is in the drying state, with (S, u) following a trajectory on ∂H dr , and the left half column (R − ) is in the wetting state, with (S, u) following a trajectory on ∂H im . These trajectories are connected at x = 0 by a horizontal jump in the (S, u) plane where u = u 0,1 ∈ (u r , u l ). This value is uniquely chosen so that the flux is continuous. The behaviour is sketched in the (S, u) plane in Fig. 5. However, it was later realized and pointed out (Raats and van Duijn 1995) that this construction fails for u l < u r . Then one has to use the scanning curves emerging from the points E l and E r as in Fig. 5. With these scanning curves, one follows the same procedure. But now the left half is in the drying state and becomes drier, while the right half is in the wetting state and becomes wetter. This is called "unconventional flow" because additionally if S l < S r then counterintuitively water flows from the dry half to the wet half. Using the Mualem model, Heinen and Raats (1999) demonstrated numerically that this type of redistribution does indeed occur.  Philip (1991) yielding "conventional" flow; (cyan) redistribution according to Raats and van Duijn (1995) yielding "unconventional" flow. The arrows indicate the direction of increasing x u S u 0,1

Scanning curves
Going one step further one can ask what happens when E i ∈ H (i = l, r ). Although redistribution results arising from interfacial area models (Pop et al. 2009) are available for general pressure and saturation initial conditions, they do not specify any directions in the induced P c -S w curves. Thus they cannot describe hysteretic redistribution in the broadest sense. This is the same reason why redistribution results for heterogeneous semi-infinite blocks (Duijn and Neef 1998) cannot be extended to cover hysteretic domains.

Remark 3
In the extended playtype model, the primary scanning curves are described by Eq. (2.3). On such curves in H one can, in principle, go back and forth, i.e. the same scanning curves are used for drying and wetting. Hence no secondary scanning curves are generated. However, if there were secondary scanning curves, they would not play a role in the horizontal redistribution cases we consider, since the right and left halves of the column can only be in one state: either drying or wetting. This applies for example to the Mualem model used by Heinen and Raats (1999).

Reformulation of the Problem
One can reduce the system of partial differential Eqs. (3.1)-(3.3) to a system of ordinary differential equations by introducing the similarity transformation: Since η → ∞ as t → 0 and x > 0, and η → −∞ as t → 0 and x < 0, the piecewise constant initial conditions (3.4), (3.5) become boundary conditions at η = ±∞. Hence substituting (4.1) into (3.1)-(3.3), using initial conditions (3.4)-(3.5) and matching conditions (3.7) one obtains the boundary value problem (P): (we will use primes to denote differentiation with respect to η) ( 4 . 4 ) (S(η), u(η)) → E r as η → ∞, (4.5) Here flux F has been redefined: it is √ t times the original flux. To obtain (4.3) we used In the formulation of Problem (P) the suction u and the flux F are continuous. Integrating by parts the first equation in (4.2), the flux continuity implies that ηS(η) is continuous as well. Hence the saturation S can only be discontinuous at η = 0. When discussing the solutions of Problem (P), we shall represent them as trajectories {(S(η), u(η)) : −∞ < η < ∞} in the (S, u) plane. The trajectories run from E l (as η → −∞) to E r (as η → ∞). In the figures below the arrows indicate the direction of increasing η (or increasing x).
Remark 4 The smoothness of the coefficients in (4.2)-(4.5) implies that the functions (S, F, u) are smooth when η = 0 and the equations are satisfied in the classical sense except at points where (S, u) moves from one of the sets H, ∂H im , ∂H dr to another.

Classification of Possible Solutions in the Half Columns R − and R +
We first consider all possible solutions in R − that satisfy boundary condition (4.4). Depending upon the location of the initial condition E l in the (S, u) plane, we distinguish three cases: Case 1. Initial condition on main wetting curve: E l ∈ ∂H im (i.e.u l = P im (S l )) Then there are four possibilities, indicated in Fig. 6. 1.4 Drying on scanning and main drying curves: There exists η 0 < 0 so that (S(η), u(η)) ∈ H and S (η) < 0 for − ∞ < η < η 0 , (S(η), u(η)) ∈ ∂H dr and S (η) < 0 for η 0 ≤ η < 0.
One may wonder why in Case 1.2 the trajectory stays on ∂H im for all η < 0. Suppose it does not. Then there exists a reversal point η * < 0 at which the trajectory switches from the main wetting curve to a drying scanning curve, as in Fig. 7 k(S(η * )) < 0, contradicting the reversal of direction at η * .
Case 3. Initial condition on main drying curve: E l ∈ ∂H dr (i.e.u l = P dr (S l )) As in Case 1 there are four possibilities, see Fig. 9.
A similar distinction of possible solutions can be made with respect to E r . We omit the details. So far we have classified all possible piecewise solutions in R − and R + . To combine the solutions in the two half columns we present one final observation.

Direction of Flow
Except for the trivial case S(η) = S l for η < 0 and S(η) = S r for η > 0 (implying that F(η) = 0 for all −∞ < η < ∞), the sign of S (η) is either strictly positive or strictly negative in the half columns R − and R + . This was demonstrated in Sect. 4.2, see also Fig. 7. We now show that this implies that the flux F cannot change its sign in the whole column: either F(η) > 0 or F(η) < 0 for all −∞ < η < ∞.
Remark 5 Redistribution and playtype hysteresis With playtype hysteresis the expression (1.5) is used. Since no time derivative of suction is involved, the redistribution problem only requires an initial saturation given by (3.4). In case of 'unconventional' flow, as suggested in Raats and van Duijn (1995), the (S, u) profile will lie on scanning curves for x = 0 which means that the saturation is constant for x < 0 and x > 0. Using this result in (3.1) we get that, F(x, t) is constant for all x ∈ R and t > 0. Since F(±∞, t) = 0 (any other value would give an unbounded suction at x = ±∞ from (3.2)) we have F(x, t) = 0 for all x ∈ R and t > 0, meaning that the suction u = u 0 = constant for all x ∈ R. This observation has the following consequence: If P im (S l ) ≤ P dr (S r ), then any horizontal connection is possible as indicated in Fig. 10. The corresponding saturation profile is 'frozen' in the sense that: for all t > 0. The suction is a undetermined constant u 0 as long as P im (S l ) ≤ u 0 ≤ P dr (S r ). If however P im (S l ) ≥ P dr (S r ) then the solution is given by classical Philip construction. Interestingly, the saturation profile remains frozen even when the domain is finite. If the domain is [−1, 1] and Neumann conditions F(±1, t) = 0 (no flow) are imposed at the boundaries for t > 0, then by the same argument one obtains saturations as in (4.7).
The cases discussed above for E l , and likewise for E r , are the building blocks in the construction of the full solution. In Sect. 4.4 we consider E l ∈ ∂H im and E r ∈ ∂H dr . In Sect. 4.5 E l and E r are arbitrarily chosen: i.e. inside H and/or on the boundaries ∂H im and ∂H dr .

Construction of Combined Solution for Initial Conditions on the Main Drying and Wetting Curves: E r ∈ ∂H dr , E l ∈ ∂H im
With reference to Fig. 11 we fix a point E r ∈ ∂H dr and let S * be such that P im (S * ) = u r . Further the curve H (S) + u = H (S r ) + u r intersects ∂H im at S = S * . Clearly 0 < S * < S * < 1. We consider five typical positions for E l with respect to the given E r .
(i) 0 < S l < S * . This is the classical Philip redistribution. In terms of Sect. 4.1 we have Case 1.2 for E l and a similar case for E r . The value of u(0) = u 0,1 , where the saturation jumps, is determined from the flux continuity. (ii) S l = S * . Special case where no flow occurs and the system is in equilibrium. Here the flux F(η) = 0 for all η ∈ R. It is Case 1.1 for E l and Case 3.1 for E r . (iii) S * < S l < S * . This is unconventional flow since the dry (left) half column becomes drier and the wet (right) half column wetter if S l < S r . Here (S(η), u(η)) ∈ H for all η ∈ R. It is Case 1.3 for E l and similar for E r . As before, the value of u(0) = u 0,3 follows from flux continuity. (iv) S l = S * . As in (iii), but now {(S(η), u(η)) : η ∈ R} belongs to a single scanning curve.
Here the saturation is continuous at η = 0. (v) S * < S l < 1. This is a rather complicated case. Depending on the position of E l with respect to E r three connections (i.e. solutions) are possible. One is shown in Fig. 11 where (S(η), u(η)) ∈ H for all η ∈ R. This situation relates to the Case 1.3 for E l and similar for E r . But there are also connections possible using part of ∂H dr (Case 1.4 for (Philip case) u l,2 = u r (Equilibrium) u l,3 < u r , S * < S l,3 < S * u l,3 < u r , S l,4 = S * u l,3 < u r , S l,5 > S * Fig. 11 Redistribution scenarios for fixed E r ∈ ∂H dr and variable E l ∈ ∂H im Fig. 12 Saturation profiles for 0 < S l < 1 when E r ∈ ∂H dr is fixed and E l ∈ ∂H im E l ) or part of ∂H im (case similar to 3.4 for E r ). These constructions will be discussed in detail in Sect. 4.5.
The saturation profiles corresponding to (i)-(v) are sketched in Fig. 12. It is relatively straightforward to check that the constructions (i)-(v) are the only ones allowed based on our discussions in Sect. 4.2.
In Pop et al. (2009) and Duijn and Neef (1998) mathematical aspects of the Philip redistribution [Case (i)] are considered. Before analysing the other cases in detail, observe that for the extended model (S(η), u(η)) ∈ H implies that H (S(η)) + u(η) is constant which means that u = − dH dS (S)S . Substituting this in (4.2) we get: Cases (iii) and (iv) are covered by (4.8) and (4.9) for η ∈ R, subject to boundary conditions S(−∞) = S l and S(+∞) = S r . Case (iv) is governed by the general theory in Duijn and Peletier (1977) whereas for Case (iii) we need to recall some of the arguments from Pop et al. (2009), Duijn andNeef (1998) and Duijn and Peletier (1977). To show how to construct the solutions in general, let us briefly consider the unconventional Case (iii). Then (S(η), u(η)) ∈ H for η ∈ R and so for the suction we have: Fig. 13 Graphical representation of (4.18). Note that since H satisfies (2.7), H (S) strictly increases with S Now consider the sub-problems u satisfies (4.10); (4.14) and u satisfies (4.11). (4.17) In these problems the saturation S − and S + will be chosen so that the suction u as well as the flux F = −D H (S)S are continuous across η = 0. From (4.10) and (4.11) it follows that (4.18) The algebraic conditions imply that if u(0) ranges from u l to u r , then S − ranges from S l to S min < S l and S + from S max > S r to S r . Boundary value problems like (P − ) and (P + ) have been studied in detail in Duijn and Peletier (1977). There it is shown that the flux at η = 0 ± depends continuously and monotonically on the boundary saturation S ± . Denoting the flux at η = 0 ± by, then F + (S + ) is continuous and strictly increasing for S r ≤ S + ≤ S max with F + (S r ) = 0, and F − (S − ) is continuous and strictly decreasing for S min ≤ S − ≤ S l with F − (S l ) = 0 (see Fig. 13).
Hence there exists a unique suction u(0) = u 0 where the fluxes intersect, as in Fig. 14. This suction uniquely determines saturations S + and S − by (4.18). Taking the composite function of solutions of (P − ) and (P + ) completes the construction of Case (iii).

Fig. 15
The redistribution scenario for general cases

Construction of Combined Solution for Arbitrary Initial Conditions: E l ,E r ∈ H ∪ ∂H im ∂H dr
For redistribution positive and negative directions are interchangeable. Therefore, without loss of generality we may assume u l < u r . (4.21) Again we fix E r , but this time E r ∈ H ∪ ∂H im ∪ ∂H dr . With assumption (4.21), we can sort out all typical cases based on the location of E l . For this purpose, and with reference to Fig. 15, the sets H 1 , H 2 and H 3 are introduced. Let S * and S * be defined as before and let the curve H (S) + u = H (P dr −1 (u r )) + u r intersect u = P im (S * ) at S = S † . The formal definitions are: Accordingly we distinguish

Case A: E l ∈ H 1
This situation is similar to Case (iii) in Sect. 4.4. Here (S(η), u(η)) ∈ H 1 ∩ H for all η ∈ R. Hence one uses Problems (P − ) and (P + ) to determine the saturations S − and S + and the . 16 a Possible switch with u l < u 0 <û; b Intersection of fluxes; c no intersection with F + (û) > F − (û) pressure u(0) for which the flux is continuous at η = 0. Figure 15 shows the trajectory running from E l = E l,1 to E r . It is comprised of part of the scanning curve through E l,1 and part of the scanning curve through E r . At η = 0 is the horizontal switch where u(0) = u 0,1 .

Case B: E l ∈ H 2
For left states in this set there are two possibilities. One is as described above for E l ∈ H 1 , that is (S(η), u(η)) ∈ H for all η ∈ R. We call this possibility 'Case B(a)'. The other possibility is more involved and needs further attention. It occurs when E l is close to ∂H dr . LetŜ denote the saturation at which the scanning curve through E l intersects ∂H dr , and letû := P dr (Ŝ) be the corresponding suction. Then forŜ ≤ S ≤ S l , (4.22) Now suppose, as in Fig. 16a, that there is a switch from the scanning curve through E l to the scanning curve through E r . Then (4.10) holds withŜ < S − < S l and in terms of the fluxes one would have a unique intersection at u(0) = u 0 satisfying u l < u 0 <û (Fig. 16b). But what if F + (û) > F − (û) as in Fig. 16c? Then the construction fails. We call this case: Case B(b).
To resolve it one needs to follow u = P dr (S) for S <Ŝ. This gives in the left column the suction saturation relation, (4.23) In the right hand column (R + ) relation (4.11) still holds. This leads to the same problem (P + ) in (R + ), but to a modified problem in R − : Fig. 17 Intersection of the fluxes withF − and F + Here the diffusivityD is in general discontinuous atŜ. But the results in Duijn and Peletier (1977) only require boundedness of the diffusivity. Thus with the notation of Sect. 4.4 (where we useF − to denote the flux in problem (P − )), we have thatF − (S − ) is strictly decreasing and continuous for S − < S l withF − (S l ) = 0. In terms of the pressure u(0), the flux is strictly increasing and continuous for u l ≤ u(0) ≤ u r withF − (u l ) = 0. Thus, as in Fig. 17, the fluxesF − and F + intersect at a unique suction u(0) = u 0 . The composite function of the solution of (P − ) and (P + ) describes the case of redistribution. Figure 18a shows the trajectory running from E l = E l,2 to E r . It is composed of part of the scanning curve through E l,2 , part of ∂H dr and then part of scanning curve through E r . The horizontal segment or switch is at η = 0 when u(0) = u 0,2 . The corresponding saturation and suction are sketched in Fig. 18b as functions of η. Due to the discontinuity in D(S) at S =Ŝ, the saturation has a kink at S =Ŝ corresponding to η = η 0 (η 0 is defined in Case 2.3, Sect. 4.2). The suction has a kink at η = 0 due to the jump in saturation.

Case C: E l ∈ H 3
There are two arrangements, denoted by Case C(a) and Case C(b), that are similar to the Cases B(a) and B(b) respectively. But in this set one has u l < P im (S * ), implying that a third type of construction is possible where part of the two scanning curves, ∂H dr as well as ∂H im are being used. We call this 'Case C(c)'.

Fig. 19 Non-intersection of the fluxes when
where the fluxes F + andF − are defined in the description of previous cases. If (4.24) holds, then the fluxes cannot intersect when u * < u(0) < u r , see Fig. 19. To resolve this situation one needs to modify Problem (P + ) by including P im , similar to the definition of Problem (P − ). The remaining argument is omitted since it is almost identical to the argument used for E l ∈ H 2 . Figure 15 shows the trajectory running from E l = E l,3 to E r .

Numerical Study
In this section we present a numerical approach for the redistribution problem and use our theoretical findings to classify and validate the computational results. For sufficiently large dimensionless column of half-length W > 0 we consider the initialboundary value problem [compare (3.1)-(3.5)] and F(±W, t) = 0 for t > 0. (5.6) For sufficiently small time t, the zero-flux boundary conditions (5.6) have a negligible influence on the redistribution process. Hence, for small t, a solution of (5.1)-(5.6) behaves as if the domain is unbounded and is close to the self-similar solutions discussed in this paper. We will reveal all typical cases.
These properties allow us to define the inverse ε = G −1 ε , so that (5.3) can be written as: Next we introduce the variable: In terms of u, v and ε , Eqs. (5.1)-(5.3) are transformed into for |x| < W and t > 0. Hence we have written (5.1)-(5.3) as a coupled system consisting of a parabolic equation for the suction u and an abstract ordinary differential equation for v. This type of splitting is well known in the mathematical literature, for instance see Schweizer (2012), van Duijn et al. (2018). Let the time interval [0, T ] be divided into N intervals of width t (T = N t) and let w n be the variable w at t = n t, with 1 ≤ n ≤ N . We calculate v n from the explicit time-discrete form of (5.11) (v 0 (x) = v(x, 0) = H (S(x, 0)) + u(x, 0)): v n = v n−1 + t ε (u n−1 , v n−1 ), for 1 ≤ n ≤ N . (5.12) We want to solve (5.10) implicitly for stability. For this, Eq. (5.10) needs to be linearized. We use a L-scheme (Pop et al. 2004) type linearization technique along with inner iterations to solve for u n . This scheme is defined by (for i = 1, 2, ....) with u 0 n = u n−1 . Here u i n is the ith iteration of the nth time step. If assumptions (A.1)-(A.4) are satisfied and if ∂ x u, ∂ u ε and ∂ v ε are bounded, then the scheme converges for L sufficiently large and for t sufficiently small (Pop et al. 2004). In (5.13) we use finite differences for spatial discretization and the whole scheme was implemented in MATLAB. In the computations we use

Numerical Results
The numerical results are obtained for k(S) = S 2 and (for simplicity) H (S) = S δ with δ = 1 40 .
For P dr and P im we took somewhat artificial expressions in order to visualize all possible cases: (1 − S) 2 and P dr (S) = 1 S − 1 + 6(1 − S) 2 .
Taking realistic (van Genuchten) expressions would make some of the cases hard to distinguish. Finally we use for ε the expression where ε = 10 −4 (fixed). Observe that ε is continuously differentiable in R. The value of γ can be chosen small as long as (S, u) ∈ H (then γ = 1), but needs a large value when (S, u) is on ∂H im or on ∂H dr . Given the size of the domain (W = 100) we show the computational results at t = 1. This is sufficiently small so that the zero-flux boundary conditions have no influence on the redistribution process. We visualize the results as trajectories in the (S, u) plane (i.e. (S(x, t), u(x, t)) where x runs from −W to W ) and as S vs. η profiles. Figure 20 shows results where E r ∈ ∂H dr is fixed and where E l ∈ ∂H im is varied. Note that the green trajectory represents unconventional flow because the trajectory moves to the left when leaving E l and when entering E r . This means that the dry half column becomes drier and the wet half column becomes wetter.  Figure 21 shows results for E l , E r ∈ H ∪ ∂H im ∪ ∂H dr . Here we took γ = 100 to ensure convergence of the iterations in (5.13). The cases discussed in Sect. 4 are accurately recovered by the computations. This validates the analysis and explains the complex behaviour of the computed saturation and suction. In fact, the agreement is excellent.

Conclusion
In this paper we discussed different hysteresis models for multiphase flow through porous media. To incorporate the effect of non-vertical scanning curves in a simple closed form, we proposed an extension to the playtype hysteresis model and showed that this model resembles the experimental scanning curves accurately. We outlined different properties of the model and discussed available physical and numerical justifications for the model.
After this we investigated horizontal redistribution in an infinite column in the context of hysteresis. In this problem, the two halves of the infinite horizontal porous column have different but constant initial saturation and pressure conditions that causes redistribution to occur. It was pointed out that existing models cannot give a complete description of the redistribution phenomenon. The extended hysteresis model was used to analyse the problem and the resulting system of equations was simplified using a similarity transformation. By distinguishing all possible cases and then using the flux and suction continuity criterion repeatedly, we constructed unique solutions for the redistribution problem. In fact, we showed that redistribution will always take place, even for unconventional cases, if the initial suction condition is different in the two halves. Moreover, we categorized all possible scenarios of redistribution into different cases.
Finally, a numerical scheme was proposed for the regularized nonlinear system of equations arising from the extended playtype hysteresis model, that converges irrespective of initial guesses. Numerical results from the scheme corroborated our analytical findings. manuscript and Prof. I. S. Pop (Hasselt University and University of Bergen) for fruitful discussions on the subject.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.