Existence and properties of solutions of extended play-type hysteresis model

This paper analyses the well-posedness and properties of the extended play-type model which was proposed in [van Duijn&Mitra (2018)] to incorporate hysteresis in unsaturated flow through porous media. The model, when regularised, reduces to a nonlinear degenerate parabolic equation coupled with an ordinary differential equation. This has an interesting mathematical structure which, to our knowledge, still remains unexplored. The existence of solutions for the non-degenerate version of the model is shown using the Rothe's method. An equivalent to maximum principle is proven for the solutions. Existence of solutions for the degenerate case is then shown assuming that the initial condition is bounded away from the degenerate points. Finally, it is shown that if the solution for the unregularised case exists, then it is contained within physically consistent bounds.


Introduction
The Richards equation is commonly used to model the flow of water and air through soil [3,21]. It is obtained by combining the mass balance equation with the Darcy's law [3,21]. In dimensionless form it reads, where the water saturation S and capillary pressure p are the primary unknowns. The saturation S is assumed to be bounded in [0, 1], k : R → R + represents the relative permeability function and g represents the gravitational acceleration assumed to be constant in our case. To close (1.1a), a relation between S and p is generally assumed. In this paper, we use the extended play-type hysteresis model (henceforth called the 'EPH' model) [16] for this purpose which equates S and p as p ∈ 1 2 (p (d) c (S) + p (i) c (S)) − 1 2 (p (d) c (S) − p (i) c (S)) sign(∂ t H(S) + ∂ t p). (1.1b) The functions p To analyse the model (1.1), we will study another class of problems, i.e.
completed with suitable boundary and initial conditions. The equivalence of (1.3) with the regularised version of (1.1) will be shown in Section 2. System (1.3) has an interesting mathematical structure as it consists of a nonlinear parabolic partial differential equation in a two-way coupling with an ordinary differential equation. This implies that standard techniques, such as the L 1 contraction [33] principle or the Kirchhoff transform [2] can not directly be applied to this system. Moreover, the maximum principle can be violated for the system (1.3) which gives rise to the overshoot phenomenon discussed in [30,44]. Further motivation for our analysis comes from establishing a physically consistent and well-posed model for hysteresis. Hysteresis refers to the phenomenon that if p is measured in an imbibition experiment where ∂ t S > 0 at all times then it follows a imbibition pressure curve denoted by p c (S)) are followed if the sign of ∂ t S changes from positive to negative or vice versa. These are referred to as the scanning curves, see Figure 1. This phenonmena was first documented in 1930 by Haines [18] and since then has been verfied by numerous experiments, [31,35,48] being some notable examples. Several different classes of models are available that include hysteresis, Lenhard-Parker model [34], interfacial area models [19,32], and TCAT (thermodynamically constrained averaging theory) models [27,27] are examples of such. Review of hysteresis from a mathematical, modelling and physical perspective can be found in [16,27,28,42].   [31]. (right) Plot from [16] showing the scanning curves for the extended play-type hysteresis model (1.1b) fitted from the left figure by tuning H(S).
However, due to the complex nature of these models, effects of hysteresis are often ignored in many practical applications and p is simply approximated by some weigthed average of p [3,21], e.g., ). The existence of weak solutions of such problems is studied in [1,2] and uniqueness is shown in [14,33] using L 1 contraction principle. However, this simplification introduces significant errors, major examples being in trapping [22], gravity fingering [37], overshoots [30] and redistribution [36]. To account for hysteresis in a comparatively simple way, the play-type hysteresis model was proposed in [5] based on thermodynamic considerations. It reads, Observe that, this gives ∂ t S ≥ 0 when p = p c (S) forces that ∂ t S = 0 which implies that the scanning curves for this model are vertical lines having constant saturation. Because of this simple structure, it can be considered to be the simplest possible model that addresses hysteresis.
The play-type hysteresis model has been studied extensively analytically due to its local and closed-form structure in contrast to other more complicated models such as the Lenhard-Parker model [34]. If the sign graph in (1.5) is regularised, then together with (1.1a), it constitutes a nonlinear pseudo-parabolic equation for S [12,25,26,41]. Existence results for such pseudo-parabolic equations can be found in [7][8][9]. The existence of solutions of the regularised play-type hysteresis model with degenerate capillary pressure and permeability functions is shown in [41]. Existence for the two-phase case is discussed in [24]. For the constant relative permeability case, the existence of weak solutions for the unregularised play-type hysteresis model is shown in [39]. In the same paper, an upscaled version of the play-type model is also derived. Uniqueness is shown in [12] for the two-phase regularised case (i.e. with dynamic capillarity). In [40] it is shown that the play-type hysteresis model does not define an L 1 -contraction. This leads to unstable planar fronts. Travelling wave solutions are investigated in [4,29,30,44]. For mathematical analysis, in many cases the pseudo-parabolic system is split into an elliptic equation coupled with an ordinary differential equation [8,9,13,25,41].
However, it was pointed out in [16] that due to approximation of the scanning curves by vertical segments, play-type hysteresis model makes certain physically inaccurate predictions. For example, it predicts that in many cases water will not redistribute when two columns having constant but different saturations are joined together. In [30] it is shown that the model predicts infinitely many interior maxima of saturation (called overshoots in this context) for high enough injection rates through a long column. However, only finite number of overshoots are observed from experiments [15]. Moreover, it is well documented that the numerical methods incorporating the play-type model become unstable if the regularisation parameter is sent to zero [10,38,44,47]. This motivated the extension of the play-type hysteresis model (EPH) given by (1.1b) in [16]. An equivalent expression was derived in [5,Eq. (35)] using thermodynamic arguments. It was used in [16] to cover all cases of horizontal redistribution and in [30] to explain the occurrence of finitely many overshoots, see also [28].
In this paper, we investigate the existence of weak solutions of the regularised EPH model and analyse the properties of its solutions. An alternative form of (1.1) is proposed in Section 2 and mathematical preliminaries are stated. In Section 3, existence of solutions for the model given by (1.3) is proven using Rothe's method. In Section 4, it is shown that a maximum principle holds for the solutions under certain assumptions. This also gives the existence of solutions for the degenerate EPH model. Section 5 is dedicated to investigating the behavior of the solutions when the regularisation parameter is passed to 0. It is shown that if a limiting solution (S, p) exists then p ∈ [p

Mathematical formulation
Let Ω ⊂ R d be a bounded open domain with ∂Ω ∈ C 1 . The L 2 (Ω) inner product in this domain is denoted by (·, ·), whereas · and · p denotes the L 2 (Ω) and L p (Ω) norms respectively for 1 ≤ p ≤ ∞. For any other space V (Ω), the norm is denoted by · V . Further, W k,p denotes the Sobolev space containing functions that have up to k th order derivatives in L p (Ω). In particular, H k (Ω) := W k,2 (Ω) and H k 0 := {u ∈ H k : u = 0 on ∂Ω in a trace sense}. Let H −1 (Ω) denote the dual of H 1 0 and ·, · the duality pairing of H 1 0 (Ω) with H −1 (Ω). The space of function having upto ℓ th order space derivatives α-Hölder continuous, will be referred to as C ℓ,α with C ℓ := C ℓ,0 .
The inequalities that are used repeatedly in our analysis include Cauchy-Schwarz inequality; Poincaré inequality, with C Ω denoting the constant appearing in it, i.e. u ≤ C Ω ∇u for u ∈ H 1 0 (Ω); Young's inequality, stating that for a, b ∈ R and σ > 0 one has and finally the discrete Gronwall's lemma, which states that if {y n }, {f n } and {g n } are non-negative sequences and y n ≤ f n + 0≤k<n g k y k , for all n ≥ 0, then y n ≤ f n + 0≤k<n f k g k exp k<j<n g j , for all n ≥ 0.
In our notation [·] + represents the positive part function, i.e.

Problem statement and assumptions
For the most part, in this paper we study the system (1.3). Completed with initial and boundary conditions, it is written as in Ω, The relation between (Ps) and EPH is established in Section 2.2. The properties of the functions D, Ψ 1/2 , u 0 and v 0 are as follows: Observe that, (2.3) combined with (A3)-(A4) imply that v(t) restricted to ∂Ω remains unchanged for all t > 0. The weak solution of (Ps) is defined as Definition 1 (Weak solution of (Ps)). The pair (u, v) with u ∈ W and v ∈ H 1 (0, T ; L 2 (Ω))∩ L 2 (0, T ; H 1 (Ω)) is a weak solution of (Ps) if u(0) = u 0 , v(0) = v 0 and it satisfies for all φ ∈ L 2 (0, T ; H 1 0 (Ω)) and ξ ∈ L 2 (Q) Due to the continuous embedding of W in C(0, T ; L 2 (Ω)), u(0) and v(0) are well-defined.
Remark 1 (Boundary conditions). For simplicity, a zero Dirichlet condition has been assumed at the boundary for our current analysis. Nevertheless, Definition 1 can be generalised to include Dirichlet, Neumann, and mixed type boundary conditions.
Remark 2 (Assumptions). The condition in (A3) that Ψ j is decreasing with respect to both the variables u and v is not required for proving the existence of the weak solutions. It is used in Section 4 to prove that the solutions are bounded in L ∞ . Similarly u 0 , v 0 ∈ L ∞ (Ω) is only used for proving the L ∞ bound. For proving the existence result in Theorem 3.1, assuming u 0 ∈ L 2 (Ω) and v 0 ∈ H 1 (Ω) is sufficient.

Relation between the regularised EPH model and (Ps)
Although (1.1) is closer to the expressions of the hysteresis models used in practice, (Ps) is more convenient to analyse mathematically. We show below that the EPH model with the sign(·) graph regularised is a particular case of (Ps). To be more precise, we start with assuming the following properties of the functions p c and k used in (1.1), (1.4) and (1.5), which are consistent with the data obtained from experiments [21,31,35].
The set of equations (1.1) cannot directly be reduced to the standard weak formulation used for partial differential equations. Thus, we consider an alternative expression to (1.1b) representing the EPH model. Completed with suitable boundary and initial conditions, the model reads Observe that, for relation (2.8) the scanning curves are given by S + b(p) = constant, instead of H(S) + p = constant as used in [16]. Moreover, if p = p The directionality imposed by hysteresis then demands that ∂ t S ≥ 0 in this case. Hence, in some open subset of Ω. Combining these observations, we assume that Here, (2.11) is the consistency criterion, a counterpart of which was stated in [16, Eq.
Remark 3 (Consistency of the scanning curves). The inequality (2.11) also guarantees that any scanning curve passing through (S 1 , p 1 ) for arbitrary S 1 ∈ (0, 1) and For the initial and boundary conditions we assume: c (S 0 ) comes from the physical constraint that p 0 stays intermediate to p  Having stated the properties of the associated functions, we now show the equivalence of the regularised (EPH) model and (Ps). For this purpose, the following transformations are introduced: (2.12) We recast (EPH) in terms of u and v. Since lim Sց0 p From these observations, we define the terms U m , V m , U M , V M that will become important later; 14) The V M > V m inequality follows from (2.13). Next, we express p = p (i) c (S) in terms of a relation between u and v. From (2.12), p = p Hence, the inverse of v i (u) exists. Let it be denoted by ρ (i) (u). In a similar way ρ (d) (u) is defined. The definitions can alternatively be summarized into From (P1) and (P3), one immediately obtains curves in the S-p plane calculated using the van Genuchten model [45]; (right) corresponding ρ (i) and .
Observe that, according to Definition 2, S has a trace on ∂Ω that does not change with time, i.e., it is fixed by S 0 .

Existence of solutions of (Pw)
The main existence result of this section is as follows: To prove this, we apply Rothe's method [23]. Let the time T be divided into N time steps of width ∆t (T = N∆t) and for any n ∈ {1, ...., N} let ∂ t w be approximated by (w n − w n−1 )/∆t for w ∈ {u, v}. Here w n stands for the value of the variable w at time t n := n∆t. The time-discrete solution is defined as Definition 3 (Time-discrete solution of (Pw)). For a given n ∈ {1, ...., N} and (u n−1 , v n−1 ) ∈ (L 2 (Ω)) 2 , the time-discrete solution of (Pw) at t = t n is a pair (u n , v n ) ∈ H 1 0 (Ω) × L 2 (Ω) which satisfies for all φ ∈ H 1 0 (Ω) and ξ ∈ L 2 (Ω), For the rest of the section the shorthand D n := D(u n , v n ), F n := F (u n , v n ) and Ψ j,n := Ψ j (u n , v n ) will be used extensively for n ∈ {1, . . . , N} and j ∈ {1, 2}. We show first that the pair (u n , v n ) exists.
From now on we assume that ∆t is small enough which guarantees the existence of solution pairs to the time discrete problems (P n ∆t ). Our goal will be to construct the solution (u, v) from the time-discrete solutions. For this purpose, we introduce the following interpolation functions: for w ∈ {u, v} and t ∈ (0, T ] the piece-wise constant interpolationsŵ ∆t ,w ∆t and the linear interpolationw ∆t are defined in Q so that for t ∈ (t n−1 , t n ] (recall that t n = n∆t), w ∆t (t) = w n ,w ∆t (t) = w n−1 ,w ∆t (t) = w n−1 + t − t n−1 ∆t (w n − w n−1 ). (3.5) As a first step we show thatū ∆t andv ∆t − v 0 are bounded uniformly in W and then we would use embedding theorems to construct the weak solution.
Proof. (Step 1) The fact that the functionsû ∆t ,ǔ ∆t ,ū ∆t belong to L 2 (0, T ; H 1 0 (Ω)) is direct. We proceed by showing their uniform boundedness. Putting φ = u n in (3.1a) and ξ = v n in (3.1b) we get, ∇u n 2 is used. Combining both inequalities and summing the results up from n = 1 to P ≤ N yields With C 0 = u 0 2 + v 0 2 and a generic constant C > 0, the discrete Gronwall's lemma is applied to get: Further, it gives two other important bounds both of which are used later, i.e.
Proof of Theorem 3.1. Lemma 3.2 shows thatū ∆t and (v ∆t − v 0 ) are bounded in W uniformly with respect to ∆t. Hence, there exists a sequence of time-steps {∆t p } p∈N with lim ∆t p = 0 such that (3.8) Due to the compact embedding of W in L 2 (Q) this further implies that, From (3.7) the following inequalities are obtained, This shows thatû ∆tp → u in L 2 (Q) as ∆t p → 0 since u −û ∆tp ≤ u −ū ∆tp + ū ∆tp − u ∆tp . By repeating this argument, one shows that the same holds forǔ ∆tp . Hence, In an identical way, for v ∆t there exists a subsequence ∆t q → 0 such that Finally, due to the bounds of ∂ tū∆t and ∂ tv∆t given in Lemma 3.2, there exists a sequence ∆t m → 0 such that, ∂ tw∆t m ⇀ ∂ t w, weakly in L 2 (0, T ; H −1 (Ω)) for w ∈ {u, v}. (3.11) We claim that (u, v) solves (Pw). Let ∆t r → 0 be a sequence that satisfies the limits (3.8)-(3.11). From (3.1) we have for φ ∈ L 2 (0, T ; H 1 0 (Ω)) and ξ ∈ L 2 (Q), Since ∂ tv∆t r ⇀ ∂ t v ∈ L 2 (Ω), the second equation directly gives (2.6b) in the limit. , v), ∇φ). Convergence of the second term (D(ǔ ∆tr ,v ∆tr )∇û ∆tr , ∇φ) remains to be shown. For this, we first observe that D(ǔ ∆tr ,v ∆tr )∇û ∆tr is bounded in L 2 (Q) uniformly with respect to ∆t. This means that a ζ ∈ L 2 (0, T ; L 2 (Ω) d ) exists such that D(ǔ ∆tr ,v ∆tr )∇û ∆tr ⇀ ζ weakly. To prove that ζ = D(u, v)∇u, we restrict the test function to φ ∈ C ∞ 0 (Q). Using the strong convergence of D(ǔ ∆tr ,v ∆tr ) and the weak convergence of ∇û ∆tr one gets with C φ = φ C 1 that as ∆t r → 0. Since the weak limit is unique, we have D(ǔ ∆tr ,v ∆tr )∇û ∆tr ⇀ ζ = D(u, v)∇u. This shows that (u, v) is a weak solution of (Pw).
4 Boundedness and existence of solutions of (EPH) 4.1 L ∞ -bounds on u and v Next, we investigate whether a solution of (Pw) satisfies the maximum principle or not. This is an interesting question primarily due to two reasons. Firstly, the maximum principle is used to prove the existence of solutions of (P EPH ) in the case when k(0) = 0. This is discussed in details later. Secondly, for pseudo-parabolic equations arising from the regularisation parameter τ , it is shown in [30,44] that the maximum principle does not hold. Having similar structure to the systems mentioned above, one might wonder if a maximum principle holds for (Pw). As it turns out, a maximum principle does hold in this case for the variables u and v under certain conditions on the advection term. This is preferred as a property of the EPH since hysteresis alone is known not to cause deviation from the maximum principle [44]. However, it is to be noted that the maximum principle does not necessarily generalise to other advective terms, as is expected in the case of pseudo-parabolic equations.
For v l := inf{v 0 } and v r := sup{v 0 } let there exist a pair {u l , u r } (u l < u r ) such that u 0 (x) ∈ [u l , u r ] for almost all x ∈ Ω and Then the weak solution (u, v) of Definition 1 satisfies v l ≤ v ≤ v r and u l ≤ u ≤ u r a.e.
The rationals behind the assumptions used in the proposition are explained in Section 4.2 in the context of (EPH).
Similarly, using the test function ξ = [v l − v] + χ t in (2.6b) yields Finally adding them yields for all t ∈ [0, T ] that

Existence of solutions of (EPH)
The main result of this section is the existence of a solution to (P EPH ). This is obtained first for the case when k(0) > 0, and then for k(0) = 0 in the absence of convective terms. (a) If k(0) > 0 then a solution of (P EPH ) in terms of Definition 2 exists.
(b) If |g| = 0, then a solution (S, p) of (P EPH ) exists for k(0) = 0. Moreover, there exists saturations 0 < S l < S r < 1 such that S(x, t) ∈ [S l , S r ] and p( To begin with, we observe that (P3) implies b ′ (p) → 0 for p → ±∞ which might cause the model to become degenerate. Therefore, we consider a non-degenerate system that approximates (P EPH ) first. Let ρ (i) (v) and ρ (d) (v), defined in (2.15), be extended to R such that For some δ > 0 small enough, let p Without loss of generality, assume that b ′ (p) > δ for p δ l < p < p δ r . Define b δ ∈ C 1 (R) to be a regularised version of b(p) such that (see Figure 3 (left)) Consequently, as p 0 satisfies (P4), p δ l < p 0 < p δ r for δ small enough, making b δ (p 0 ) = b(p 0 ). The values v l , v r , u l , u r , S l and S r are shown, as well as the set (v 0 , u 0 ).
Proof of Theorem 4.1(a). Define S δ := v δ −u δ . From Lemma 4.1 it follows that there exists a sequence {δ r } r∈N with lim δ r = 0 such that (→ implies strong and ⇀ implies weak convergence) Hence, following the proof of Theorem 3.1 we conclude that u, v, p and S satisfy for all φ ∈ L 2 (0, T ; H 1 0 (Ω)) and ξ ∈ L 2 (Q). For completeness, we still need to show that u = b(p) and S ∈ [0, 1] a.e. To show the former, observe that b(p δr ) → u in L 2 (Q), since, 3) The first term on the right vanishes as u δ → u strongly. For the second we use (Pδ), giving b δr (p δr ) − b(p δr ) ≤ Cδ r [1 + p δr ] for some C > 0, which approaches 0 as δ r → 0. Now, from (P3) one gets As δ r → 0, the first term in the right vanishes due to the weak convergence of p δr and the second term vanishes due to (4.3), thus proving u = b(p).
To see this, we use the test function ξ = [V m − v] + in (4.2b) to get Here, as v < V m and u < U Similarly one shows that S ≤ 1 which concludes our proof.
5 Behaviour when τ → 0 In this section we address one other important question: what is the behavior of (S, p) for the limit τ → 0 which corresponds to the unregularised EPH model. Ideally p ) should be satisfied in the pure hysteresis case, i.e., when τ → 0. One might also wonder whether a limiting solution exists in this case or not. The problem is addressed for the play-type hysteresis model (1.5) in [39, Theorem 3.2] for a reduced case (linear, no advection etc.), however, with stochastic variation of coefficients included. For the nonlinear case, no results are available to our knowledge even for the play-type model. In this section, we take a step towards understanding the limiting case τ → 0. Proposition 2. Let the assumptions of Theorem 4.1 hold. For a given τ > 0, let (S τ , p τ ) denote a weak solution of (P EPH ) in terms of Definition 2. Further let u τ := b(p τ ) and v τ := S τ + b(p τ ). Then, for a constant C > 0 independent of τ , one has Proof. Observe that Φ τ can be rewritten simply as We consider here the case when |g| = 0 and k(0) = 0, the case |g| = 0 and k(0) > 0 being simpler. From Theorem 4.1(b), S τ satisfies 0 < S l ≤ S τ ≤ S r < 1 a.e. in Q, where S l , S r are independent of τ . For D(u, v), defined as in (2.19), this implies that there exists D m = k(S l )/ max p∈R {b ′ (p)} > 0 such that (A1) is satisfied. The pair (u τ , v τ ) satisfies This implies that if a limit (u, v) exists of (u τ , v τ ) as τ → 0, then u ∈ [ρ (i) (v), ρ (d) (v)] in Q implying that the pressure in this limit stays within the bounds of the hysteretic pressure curves. Another immediate consequence of Proposition 2 is the boundedness of T 0 ∇u τ 2 as evident from (5.5). This implies Corollary 1. Under assumptions of Proposition 2, there exists u ∈ L 2 (0, T ; H 1 0 (Ω)) ∩ L ∞ (0, T ; L 2 (Ω)) such that u τ ⇀ u in L 2 (0, T ; H 1 0 (Ω)) weakly as τ → 0.

Conclusion
In this paper, we analysed the extended play-type hysteresis (EPH) model, proposed in [16], for the unsaturated flow case. A modification of the model is proposed that conforms with the standard weak formulation used in analysis. A regularised version of the problem is considered and it is shown to be equivalent to a nonlinear parabolic equation coupled with an ordinary differential equation. Using Rothe's method, the existence of weak solutions is first proven for the equivalent system and then for the (EPH) model. Solutions are shown to exist even when the capillary pressure functions blow up and relative permeability vanishes for zero saturation. To cover the latter case, the maximum principle needs to hold, which is proven to be the case in the absence of an advection term. The consequences of passing the regularisation parameter to zero are explored, which show that the pressure stays bounded by capillary pressure curves in the limit.