Zeitschrift für Analysis und ihre Anwendungen

Journal for Analysis and its Applications Volume 19 (2000), No. 4, 977–997

Hysteresis in Filtration through Porous Media

F. Bagagiolo and A. Visintin

Abstract. We study an evolution problem for filtration through porous media, accounting for hysteresis in the saturation versus pressure constitutive relation. We provide a weak formulation of the problem, assuming that the memory effect in the constitutive relation consists not only of a rate-independent component but also of a rate-dependent one. We prove an existence result, which also applies to the case where the hysteresis operator is of Preisach-type.

Keywords: Degenerate parabolic equations, hysteresis, filtration, variational inequalities AMS subject classification: Primary 35K65, 46J50, 76S05, secondary 49J40

1. Introduction

Although a large technical literature accounts for hysteresis effects in porous media filtration, apparently the analytical aspects of this phenomenon have not yet been studied. A simplified model leads to the system

\begin{matrix}φ\frac{∂s}{∂t} − ∇ · k(∇u + ρgz⃗) = 0  in D′(Q) \\ s = F̃(u), k = k̃(s)  in Q \end{matrix} \} .  (1.1)

This must be coupled with appropriate initial and boundary conditions, including a seepage condition of Signorini type. The saturation s and the pressure u are unknown. The quantity

φ ∈ [0, 1] represents the porosity k the hydraulic conductivity g the gravity acceleration ρ the density of the fluid z⃗ the upward vertical unit vector.

The function k̃ : [0, 1] → ℝ^{+} is prescribed and continuous.

Equation (1.1)_{1} follows from the mass balance and Darcy’s law. The dependence of s upon u is formally represented by the operator F̃; the description of the latter is one of the main issues of this work. Experimental evidence indicates the occurrence of a quantitatively relevant hysteresis effect, which has occasionally been represented by Preisach-type models in the engineering literature.

F. Bagagiolo: Univ. of Padova, Dept. Pure & Appl. Math., Via Belzoni 7, 35131 Padova, Italy A. Visintin: Univ. of Trento, Dept. Math., Via Sommarive 14, 38050 Povo-Trento, Italy ISSN 0232-2064 / $ 2.50 \begin{matrix}c \\ ◯ \end{matrix} Heldermann Verlag Berlin

This problem exhibits several interesting features:

(i) The parabolic part of equation (1.1)_{1} may degenerate, as s is bounded. (ii) The elliptic part of equation (1.1)_{1} may also degenerate, as k may vanish. (iii) The s versus u constitutive relation exhibits hysteresis. (iv) The coefficient k depends on a hysteresis-dependent term, s.

In this paper we discuss the hysteresis relation, and propose a model for which we are able to prove existence of a weak solution. In Section 2 we briefly illustrate the equations representing saturated-unsaturated filtration through porous media. In Sections 3 and 4 we then discuss the s versus u constitutive relation, propose an equation which combines hysteresis and time relaxation, and study its analytical properties. In Section 5 we formulate an initial- and boundary-value problem in the framework of Sobolev spaces for system (1.1), and prove existence of a weak solution. Finally, in Section 6 we briefly discuss the Preisach model, which is used to represent the hysteretic component of the s versus u relation.

We refer to Bear [5] and to Fredlund and Rahardjo [8] for a presentation of the physical and engineering background; we also refer to Poulovassilis and Childs [18], Poulovassilis and Tzimas [20], Poulovassilis and El-Ghamry [19], Mualem [15, 16], Kool and Parker [11] and to the references therein for experimental studies of the hysteresis relation. Saturated-unsaturated flow in porous media with free boundary has been studied in a large number of papers, in particular we refer to works of Baiocchi [4], Torelli [21], Alt [1], and Gilardi [9]. A model analogous to (1.1), with a fairly general saturation versus pressure constitutive relation but without hysteresis, has been studied by Alt, Luckhaus and Visintin [2], Otto [17] and by Bagagiolo [3]. We refer to the monographs of Krasnosel’skiĭ and Pokrovskiĭ [12], Mayergoyz [14], Visintin [22], Brokate and Sprekels [6] and Krejčí [13] for hysteresis models.

2. Filtration in a porous medium

Let Ω be a bounded, open and connected set of ℝ^{3}, with Lipschitz boundary, representing the region occupied by the porous medium (see Fig. 1). Let the boundary of

Figure 1: Section of a porous dam with two reservoirs Ω be divided in two parts, namely Γ_{1} the impermeable part and Γ_{2} the part in contact with either water or air. We assume that Γ_{1} and Γ_{2} are Lipschitz bidimensional manifolds. Let [0, T] be a time interval, with T > 0, and define Q = Ω × (0, T), Σ_{1} = Γ_{1} × (0, T) and Σ_{2} = Γ_{2} × (0, T). Let us denote by s, u and k the saturation of the medium, the pressure of water inside the medium and the hydraulic conductivity of the medium, respectively.

We suppose that we are in the range of validity of Darcy’s law, which is essentially an experimental law. That law yields the following relation between the flux q of water inside the porous medium, pressure and hydraulic conductivity

q = −k∇(u + ρgz)  in Q,  (2.1)

where z is the vertical coordinate of the point x, g is the gravity acceleration and ρ is the density of the fluid (i.e. water). For the sake of simplicity, let us omit the porosity coefficient φ of the medium. From (2.1) and using the equation of continuity for the content of water inside any closed region of the medium, we obtain the equation

s_{t} − ∇ · [k∇(u + ρgz)] = 0  in Q,  (2.2)

where s_{t} is the time derivative of s, ∇ is the spatial gradient, and “∇·” is the divergence operator.

Let P be a non-negative function defined on Σ_{2}, representing the datum for the pressure u. Typically P vanishes on the part of Σ_{2} in contact with air, whereas it coincides with the corresponding hydrostatic pressure of the reservoir on the part of Σ_{2} in contact with water. Let us denote by ν the outward normal unit vector to Ω. On Σ_{1} ∪ Σ_{2} we have the following conditions:

k∇(u + ρgz) · ν =0  on Σ_{1}  (2.3)
u^{+} =P  on Σ_{2}  (2.4)
k∇(u + ρgz) · ν ≤0  on Σ_{2} ∩ {u = 0}  (2.5)
k∇(u + ρgz) · ν =0  on Σ_{2} ∩ {u < 0}.  (2.6)

By (2.1), (2.3) means that there is no flux through the impervious part Σ_{1}; (2.4) means that the positive part of the pressure is prescribed on Σ_{2}; (2.5) means that through the part of Σ_{2} where the pressure vanishes, that is where the medium is in contact with air, water can only flow outward; (2.6) means that through the part of Σ_{2} where the pressure is negative (the boundary of the so-called capillary fringe) there is no flux.

Conditions (2.5) and (2.6), together with (2.4), are equivalent to the variational inequality of “Signorini type”

k∇(u + ρgz) · ν(u − v) ≤ 0 on Σ_{2}, ∀v : Σ_{2} → ℝ such that v^{+} = P.  (2.7)

We must also prescribe an initial condition for the saturation s

s(·, 0) = s^{0}(·)  in Ω.  (2.8)

The constitutive relation between the saturation s and the pressure u is typically represented by a relation of the form

s(x, t) ∈ f(u(x, t)),  ∀ (x, t) ∈ Q,  (2.9)

where f : ℝ → [0, 1] is a maximal monotone graph as in Figure 2.

Figure 2: Non-hysteretic saturation versus pressure constitutive relation

The value ŝ represents a possible irreducible level of saturation. Moreover, the hydraulic conductivity is represented by a nonnegative continuous function k(s) of the saturation, as in Figure 3; in particular k > 0 in (ŝ, 1].

Figure 3: Hydraulic conductivity versus saturation constitutive relation

In the sequel we shall discuss and amend the constitutive relation (2.9).

3. The saturation versus pressure constitutive relation

Experimental evidence (see Figure 4) indicates that, at any point x ∈ Ω, s(x, t) depends not only on u(x, t), but also on the initial value s^{0}(x) (see also Remark 6.5) and on the previous evolution of u at the same point, u(x, ·). We assume that here x occurs just as a parameter, and for the moment we do not display it in the study of the constitutive behavior.

At first, let us then consider a dependence of the form

s(t) = [F̃(u, s^{0})](t)  ∀ t ∈ [0, T].  (3.1)

Here u represents the whole function u : [0, T] → ℝ. Although s(t) cannot depend on u_{|(t,T]}, for formal reasons it is convenient to deal with the whole function, and then explicitly require that at time t, [F̃(u, s^{0})](t) does not depend on u_{|(t,T]} (causality). The operator F̃ is then called a memory (or Volterra) operator.

Figure 4: Hysteretic constitutive relation

On account of experimental evidence, as a first hypothesis it seems natural to assume that F̃ is a hysteresis operator. By this we mean that F̃ is rate-independent, that is, for any non-decreasing homeomorphism φ : [0, T] → [0, T],

F̃(u ∘ φ, s^{0}) = F̃(u, s^{0}) ∘ φ  in [0, T].  (3.2)

According to the current terminology, this characterizes purely hysteretic effects. We can assume that F̃ is continuous from C^{0}([0, T]) × ℝ to C^{0}([0, T]), consistently with most of the known models of hysteresis.

When coupling this constitutive relation with the partial differential equation (1.1)_{1} of the Introduction, we must insert the dependence on the parameter x ∈ Ω:

s(x, t) = [F̃(u(x, ·), s^{0}(x))](t)  ∀ (x, t) ∈ Q.  (3.3)

Difficulties arise in proving existence of a solution for the corresponding initial- and boundary-value problem. By a standard procedure we might approximate the problem, derive a priori estimates, then try to pass to the limit. On account of the occurrence of the memory operator, it seems especially convenient to use implicit time discretization. Let us denote the approximation parameter by m ∈ ℕ, and the approximate solution by (u_{m}, s_{m}). Uniform estimates for u_{m} in L^{2}(0, T; H^{1}(Ω)) can be derived by multiplying the approximate equation by u_{m}, and then integrating over Q. But this does not yield convergence of u in C^{0}([0, T]) a.e. in Ω (not even for a subsequence), hence it does not suffice to pass to the limit in the memory operator. In order to derive stronger a priori estimates, we might try to multiply the approximate equation by \frac{∂u_{m}}{∂t}, and then integrate over Q. However, difficulties arise in dealing with the elliptic term.

If the dependence of s on u were without memory, s = f(u) with f ∈ C^{0}(ℝ) say, then we would apply the Kirchhoff transformation:

K : u ↦ ũ := ∫_{0}^{u}(k̃ ∘ f)(ξ) dξ,  (3.4)

so that ∇ũ = k̃(f(u))∇u a.e. in Q. As k̃ and f are non-decreasing, K would be invertible; (3.4) would then be equivalent to ũ ∈ α(u), with α a (possibly multi-valued) maximal monotone operator. We might then write the elliptic part as −Δũ, and express (1.1)_{1} in terms of s and ũ. This procedure was used in Alt, Luckhaus and Visintin [2] to deal with the problem without hysteresis. However, it is clear that the Kirchhoff transformation cannot be applied whenever memory occurs in the s versus u constitutive relation.

These difficulties induce us to revise the formulation of the model. Although we are not able to derive a uniform estimate on the pressure rate, we conjecture that this rate should not be too large, even on the (rather slow) time scale typical of filtration phenomena. We then propose to insert in the s versus u constitutive relation a term which penalizes high rates. By this we shall account not only for hysteresis but also for a small rate-dependent component of the memory.

Let us detail the construction of our model. We suppose that the hysteresis branches (hysteresis loops) occur only for values of u which belong to a bounded set, say [u_{1}, u_{2}] ⊂ ℝ; that is on (−∞, u_{1}] ∪ [u_{2}, +∞] the operator F̃ acts just as a superposition operator: s ∈ ϕ̃(u) where ϕ̃ is a (possibly multi-valued) maximal monotone graph, with lim_{u→−∞} ϕ̃(u) = ŝ and lim_{u→+∞} ϕ̃(u) = 1 (see Figure 4). Hence we represent the hysteresis relation by an operator

F : C^{0}([0, T]; [u_{1}, u_{2}]) × [s_{1}, s_{2}] → C^{0}([0, T]; [s_{1}, s_{2}]),

where [s_{1}, s_{2}] ⊆ [ŝ, 1]. Let τ and σ be the truncation functions

τ(ξ) = \{ \begin{matrix}u_{1} & if ξ ≤ u_{1} \\ ξ & if u_{1} ≤ ξ ≤ u_{2} \\ u_{2} & if ξ ≥ u_{2} \end{matrix}  and  σ(ξ) = \{ \begin{matrix}s_{1} & if ξ ≤ s_{1} \\ ξ & if s_{1} ≤ ξ ≤ s_{2} \\ s_{2} & if ξ ≥ s_{2}. \end{matrix}

For any s^{0} ∈ [ŝ, 1], we define the hysteresis operator F(τ(·), σ(s^{0})) which acts on a continuous function u : [0, T] → ℝ in the following way. For any t ∈ [0, T] such that u(t) ∈ (u_{1}, u_{2}), let us define

t̃= min {t′ ∈ [0, t] : u(t′′) ∈ [u_{1}, u_{2}] ∀ t′′ ∈ [t′, t]}.

For any t ∈ [0, T], we then have [F(τ(u), σ(s^{0}))](t)

= \{ \begin{matrix}s_{1} & if u(t) ≤ u_{1} \\ [F(u(· + t̃), s_{1})](t − t̃) & if u(t) ∈ (u_{1}, u_{2}), t̃> 0, u(t̃) = u_{1} \\ [F(u(· + t̃), s_{2})](t − t̃) & if u(t) ∈ (u_{1}, u_{2}), t̃> 0, u(t̃) = u_{2} \\ [F(u, σ(s^{0}))](t) & if u(t) ∈ (u_{1}, u_{2}), t̃= 0 \\ s_{2} & if u(t) ≥ u_{2} \end{matrix} (3.5)

where u(· + t̃) : [0, t − t̃] → [u_{1}, u_{2}] and t′ ↦ u(t′ + t̃).

Hence, for a suitable maximal monotone graph

ϕ : ℝ → [ŝ − s_{1}, 1 − s_{2}],

satisfying ϕ(ξ) = 0 for every ξ ∈ [u_{1}, u_{2}] and

\begin{aligned}lim_{ξ→−∞} ϕ(ξ) & = ŝ − s_{1} \\ lim_{ξ→+∞} ϕ(ξ) & = 1 − s_{2} \end{aligned}

the constitutive relation (3.1) can be written as

s(t) ∈ [F̃(u, s^{0})](t) := [F(τ(u), σ(s^{0}))](t) + ϕ(u(t)).  (3.6)

Under a natural assumption, so-called piecewise monotonicity (namely, the monotonicity of all the hysteresis branches, see Visintin [22: p. 62]), for any choice of the initial state s^{0} ∈ [s_{1}, s_{2}], we can invert the operator F(·, s^{0}). If F(·, s^{0}) is a hysteresis operator, then the same holds for its inverse G(·, s^{0}).

Now let us define the following aximal monotone graph on [ŝ, 1]:

β(s) = \{ \begin{matrix}ϕ^{−1}(ŝ − s) − u_{1} & if s ≤ s_{1} \\ 0 & if s_{1} ≤ s ≤ s_{2} \\ ϕ^{−1}(1 − s) − u_{2} & if s ≥ s_{2}. \end{matrix}

We define the operator

[G̃(s, s^{0})](t) = [G(σ(s), σ(s^{0}))](t) + β(s(t)),

where G(σ(·), σ(s^{0})) is defined in a similar way as F(τ(·), σ(s^{0})). It is easy to check that G̃ = F̃^{−1}. Hence (3.6) can be written in the equivalent form (see Figure 5)

u(t) ∈ [G̃(s, s^{0})](t) = [G(σ(s), σ(s^{0}))](t) + β(s(t))  (3.7)

We then insert a rate-dependent memory effect, and write our constitutive relation in the form

u(t) ∈ [G̃(s, s^{0})](t) + α\frac{ds}{dt}  (3.8)

where α is a small positive relaxation constant.

In the next section we show that under natural assumptions this equation defines a continuous operator F̃_{α} : L^{2}(0, T) → H^{1}(0, T) : u ↦ s, for any fixed s^{0} ∈ [ŝ, 1]. Obviously, F̃_{α} is rate-dependent. This relation is then extended to the space-distributed problem by inserting the dependence on the parameter x.

The regularity properties of the operator F̃_{α} will allow us to prove existence of a solution of the modified problem. It would then be natural to consider the behavior of the solution of our problem as α → 0. But, as it might be expected, in this limit we encounter the same difficulties that we pointed out above for the purely hysteretic constitutive relation.

4. Analytical properties of the relaxed constitutive relation

As in the previous section, we represent the hysteresis operator G̃ (see (3.7)) as the sum of two operators. One of them, denoted by G, is bounded and accounts for hysteresis; the other one, denoted by β, is an unbounded maximal monotone graph on ℝ with domain [ŝ, 1] (see Figure 5). The graph β may contain horizontal segments, which may correspond to vertical segments (jumps) of the s versus u relation (if α = 0 this is a free boundary problem).

Figure 5: Hysteretic and monotone parts of the (inverse) constitutive relation

For the moment, let the time t ∈ [0, T] be the only meaningful independent variable, and let us regard x just as a parameter. Accordingly here we drop the notation x among the entries of all functions.

We recall the reader that, for any fixed s^{0} ∈ ℝ,

G(·, s^{0}) : C^{0}([0, T]) → C^{0}([0, T])

is called a hysteresis operator whenever it is causal and rate-independent (see the previous section and also Section 6). We consider the following ordinary differential inclusion:

\begin{matrix}α\frac{ds}{dt} + G(s, s^{0}) + β(s) ∋ u  a.e. in [0, T] \\ s(0) = s^{0} \end{matrix}  (4.1)

where s^{0} ∈ [ŝ, 1] (with 0 < ŝ < 1) is fixed, α is a positive real number and u is a given function.

Proposition 4.1. Let β be an unbounded maximal monotone graph on ℝ with domain [ŝ, 1], and G(·, s^{0}) : C^{0}([0, T]) → C^{0}([0, T]) be a bounded and Lipschitz continuous hysteresis operator (in particular a causal operator). Then, for any u ∈ L^{2}(0, T), there exists a unique solution s ∈ H^{1}(0, T) of (4.1). This defines the solution operator

F : L^{2}(0, T) → H^{1}(0, T) : u ↦ s.  (4.2)

Proof. Let us take u ∈ L^{2}(0, T). To solve (4.1) we use an implicit time discretization. For every m, n ∈ ℕ \ \{0\}, n ≤ m, let us define h = \frac{T}{m}, u_{m}^{n} = \frac{1}{h} ∫_{(n−1)h}^{nh} u(t) dt and s_{m}^{0} = s^{0}. Let us assume that, for every 0 ≤ l < n, the values s_{m}^{l−1} are known. By the memory properties of G, for every continuous function s which is linear on the intervals [lh, (l + 1)h] and such that s(lh) = s_{m}^{l} for every l, the value of [G(s, s^{0})]((l + 1)h) depends only on s((l + 1)h). Hence, for every ξ ∈ ℝ, we can define ψ_{m}^{n}(ξ) as the value of [G(s, s^{0})]((l + 1)h) where s is a function as above such that s((l + 1)h) = ξ. It is easy to prove that the function ψ_{m}^{n} is bounded and Lipschitz continuous on ℝ, uniformly with respect to m and n, as the same properties hold for G.

Let us suppose that s_{m}^{n−1} ∈ [ŝ, 1] is known, and consider the following problem: to find s_{m}^{n} ∈ [ŝ, 1] such that

α\frac{s_{m}^{n} − s_{m}^{n−1}}{h} + ψ_{m}^{n}(s_{m}^{n}) + β(s_{m}^{n}) ∋ u_{m}^{n}.  (4.3)

We use a fixed point procedure to solve (4.3). Let us take r ∈ [ŝ, 1] and consider the problem: to find ξ_{r} ∈ [ŝ, 1] such that

α\frac{ξ_{r} − s_{m}^{n−1}}{h} + ψ_{m}^{n}(r) + β(ξ_{r}) ∋ u_{m}^{n}.  (4.4)

Let b be a convex and lower semicontinuous primitive of β, that is ∂b = β (the subdifferential of β); in particular, b ≡ +∞ in (−∞, ŝ[∪]1, +∞). Problem (4.4) can be solved by minimizing on ℝ the continuous strictly convex function

J(ξ) = \frac{α}{2h}ξ^{2} + (ψ_{m}^{n}(r) − \frac{αs_{m}^{n−1}}{h} − u_{m}^{n})ξ + b(ξ),  (4.5)

and using the fact that ξ_{r} minimizes J if and only if 0 ∈ ∂J(ξ_{r}).

By the coercivity and strict convexity of J, the minimizer is unique and belongs to [ŝ, 1]. Hence we get a function r ↦ ξ_{r} which maps [ŝ, 1] onto itself. We claim that J is a contraction if h is sufficiently small; hence it has a unique fixed point, which is the unique solution of problem (4.3). In fact, let us take r, s ∈ [ŝ, 1] and suppose that ξ_{r} < ξ_{s}. By (4.4) and using the monotonicity of β, we obtain

0 < ξ_{s} − ξ_{r} ≤ \frac{h}{α}(ψ_{m}^{n}(r) − ψ_{m}^{n}(s))

from which, recalling that ψ_{m}^{n} is Lipschitz, the conclusion easily follows.

Let us now denote by s_{m} the linear time interpolate of the nodal values s_{m}^{n}. Obviously, s_{m} ∈ L^{∞}(0, T) and ∥s_{m}∥_{L^{∞}(0,T)} ≤ 1 for every m. If we multiply (4.3) by s_{m}^{n} − s_{m}^{n−1} and sum over n then, using the fact that ∂b = β, for every 0 < l ≤ m we get

hα ∑_{n=1}^{ℓ} (\frac{s_{m}^{n} − s_{m}^{n−1}}{h} )^{2} + h ∑_{n=1}^{ℓ} \frac{s_{m}^{n} − s_{m}^{n−1}}{h} [ψ_{m}^{n}(s_{m}^{n}) − u_{m}^{n}] ≤ b(s^{0}) − b(s_{m}^{ℓ}) ≤ C (4.6)

where the last inequality, independent on ℓ and m, holds since −b is upper semicontinuous and s_{m}^{ℓ} ∈ [ŝ, 1]. Let us denote by \overline{ψ}_{m} (respectively, \overline{u}_{m}, \overline{s}_{m}) the piecewise constant interpolate of the nodal values ψ_{m}^{n}(s_{m}^{n}) = [G(s_{m}, s^{0})](nh) (respectively u_{m}^{n}, s_{m}^{n}). By the properties of G, we conclude that \overline{ψ}_{m} → G(s) strongly in L^{2}(0, T), moreover we have also the strong convergence \overline{u}_{m} → u in L^{2}(0, T).

We can write (4.6) in the form

α ∫_{0}^{ℓh} [(\frac{ds_{m}}{dt} )^{2} + \frac{ds_{m}}{dt} (\overline{ψ}_{m} − \overline{u}_{m})]dt ≤ C.  (4.7)

By (4.7) we get ∥s_{m}∥_{H^{1}(0,T)} ≤ C for all m ∈ ℕ \ {0}. Hence there exists s ∈ H^{1}(0, T) such that s_{m} → s weakly in H^{1}(0, T) whence, by compactness, s_{m} → s strongly in C^{0}([0, T]). Moreover, s(t) ∈ [ŝ, 1] for all t and s(0) = s^{0}.

Now we show that s is the solution of problem (4.1). Note that (4.1) is equivalent to

(u(t) − α\frac{ds}{dt}(t) − G(s, s^{0})(t))(s(t) − ξ) ≥ b(s(t)) − b(ξ)  (4.8)

for all ξ ∈ ℝ, a.e. in [0, T]. Let us take ξ ∈ ℝ and φ ∈ C^{0}([0, T]), φ ≥ 0. From (4.3) we get

∫_{0}^{T} φ(t)(\overline{u}_{m} − α\frac{ds_{m}}{dt} − \overline{ψ}_{m})(\overline{s}_{m} − ξ) dt ≥ ∫_{0}^{T} φ(t)(b(\overline{s}_{m}) − b(ξ)) dt.  (4.9)

Passing to the inferior limit m → +∞ in (4.9), using the lower semicontinuity of b, by the arbitrariness of φ we obtain (4.8) and hence (4.1).

Finally, let us prove the uniqueness. Let s_{1} and s_{2} be solutions of problem (4.1). We multiply the difference of the two inclusion by s_{1} −s_{2} and integrate over [0, t]. Using the monotonicity of β and the Lipschitz continuity of G, we obtain

\frac{α}{2}|(s_{1} − s_{2})(t)|^{2} ≤ L∥s_{1} − s_{2}∥_{C^{0}([0,t])} ∫_{0}^{t} |(s_{1} − s_{2})(τ)| dτ  (4.10)

where L is the Lipschitz constant of G. From (4.10), the inequality

\frac{α}{2}∥s_{1} − s_{2}∥_{C^{0}([0,t])} ≤ L ∫_{0}^{t} ∥s_{1} − s_{2}∥_{C^{0}([0,τ])}dτ

easily follows and hence, using the Gronwall lemma, we get s_{1} = s_{2} in [0, T]

Proposition 4.2. The solution operator F, defined in (4.2), is Lipschitz continuous from L^{2}(0, T) to C^{0}([0, T]).

Proof. Let s_{1} and s_{2} be two solutions of problem (4.1) with different sources u_{1} and u_{2}, respectively. Let us define s̃ = s_{1} − s_{2} and ũ = u_{1} − u_{2}. Subtracting one inclusion from the other, multiplying by s̃ and integrating, we get, in a similar way as before,

\frac{α}{2}|s̃(t)|^{2} ≤ L∥s̃∥_{C^{0}([0,t])} ∫_{0}^{t} |s̃(τ)| dτ + ∫_{0}^{t} |ũ(τ)| |s̃(τ)| dτ  ∀ t ∈ [0, T].  (4.11)

Note that, since the right-hand side of (4.11) is monotone in t, the same inequality holds even if we replace, in the left-hand side only, the instant t with any other instant \overline{t} ∈ [0, t] and hence if we replace the left-hand side with \frac{α}{2}∥s̃∥_{C^{0}([0,t])}^{2}. Moreover, for every ε > 0 we can majorize the two terms in the right-hand side by

\begin{matrix}εL^{2}T ∥s̃∥_{C^{0}([0,T])}^{2} + \frac{1}{ε} ∫_{0}^{t} |s̃(τ)|^{2}dτ \\ \frac{1}{ε} ∫_{0}^{t} |ũ(τ)|^{2}dτ + ε ∫_{0}^{t} |s̃(τ)|^{2}dτ, \end{matrix}

respectively. Hence, if ε > 0 is suitably small, we obtain

0 ≤ (\frac{α}{2} − T (εL^{2} + ε))∥s̃∥_{C^{0}([0,t])}^{2} ≤ \frac{1}{ε} ∫_{0}^{t} ∥s̃∥_{C^{0}([0,τ])}^{2}dτ + \frac{1}{ε}∥ũ∥_{L^{2}(0,T)}^{2}.  (4.12)

By (4.12), using the Gronwall inequality, we easily get the Lipschitz continuity of the operator F

Inclusion of the spatial dependence.

Now we insert into (4.1) the dependence on x ∈ Ω. First of all, let us consider a hysteresis operator

G̃: C^{0}([0, T]) × ℝ → C^{0}([0, T]).

For every function w ∈ L^{2}(Ω; C^{0}([0, T])) and any initial state ξ^{0} ∈ L^{2}(Ω), we define

[G(w, ξ^{0})](x, t) = [G̃(w(x, ·), ξ^{0}(x))](t)  a.e. in Ω, ∀ t ∈ [0, T].  (4.13)

Roughly speaking, we apply the operator G̃ at almost every x ∈ Ω.

It is easy to see that G is causal, that is if w_{1} = w_{2} in [0, t] a.e. in Ω and ξ_{1}^{0} = ξ_{2}^{0} a.e. in Ω, then

[G(w_{1}, ξ_{1}^{0})](·, t) = [G(w_{2}, ξ_{2})](·, t) a.e. in Ω.

Moreover, from the continuity properties of G̃, it follows that (see Visintin [22: p. 258]) [G(w, ξ^{0})](x, ·) ∈ C^{0}([0, T]) a.e. x ∈ Ω and that G(w, ξ^{0}) : Ω → C^{0}([0, T]) is measurable and satisfies the following property: if w_{n} → w uniformly in [0, T] a.e. in Ω and ξ_{n}^{0} → ξ^{0} a.e. in Ω, then G(w_{n}, ξ_{n}^{0}) → G(w, ξ^{0}) uniformly in [0, T] a.e. in Ω. Hence, we can reformulate problem (4.1) as

α\frac{∂s}{∂t}(x, t) + [G(s, s^{0})](x, t) + β(s(x, t)) ∋ u  a.e. in (x, t) ∈ Q;  (4.14)

here u ∈ L^{2}(0, T; H^{1}(Ω)) and s^{0} ∈ L^{∞}(Ω) are assumed to be given.

In order to find a solution of problem (4.14), we solve (4.1) for almost every x ∈ Ω. Let us call F_{m}^{n} the function on ℝ such that, for every ξ ∈ ℝ, F_{m}^{n}(ξ) is the unique solution of the discrete problem (4.3), where u_{m}^{n} is replaced by ξ. We claim that for sufficiently small h, F_{m}^{n} is Lipschitz continuous and its Lipschitz constant is independent on s_{m}^{n−1} ∈ ℝ. In fact, if s = F_{m}^{n}(ξ) and w = F_{m}^{n}(ζ), we subtract the two respective inclusions (4.3) and using the monotonicity of β and the Lipschitz continuity of ψ_{m}^{n}, we get the conclusion by a procedure similar to the one we used in the proof of Proposition 4.1. Hence, if u ∈ L^{2}(0, T; H^{1}(Ω)), inserting the parameter x ∈ Ω, for every m and n we obtain s_{m}^{n} ∈ H^{1}(Ω).

5. Weak formulation of the problem and existence of a solution

Let Ω ∈ ℝ^{3} be open, bounded and connected, and Γ ⊆ ∂Ω a Lipschitz manifold, with positive bidimensional measure. Let [0, T] be a time interval and let us set Q = Ω×(0, T) and Σ = Γ × (0, T).

We consider the functions P ∈ C^{0}([0, T], H^{1}(Ω)) ∩ H^{1}(0, T; L^{2}(Ω)), P ≥ 0 a.e. in Q, and s^{0} ∈ L^{∞}(Ω). We define the convex set

K = {u ∈ L^{2}(0, T; H^{1}(Ω)) : (γ_{0}u)^{+} = P on Σ}

where γ_{0} is the trace operator: H^{1}(Ω) → H^{\frac{1}{2}} (Γ). Let G be a hysteresis operator and β be a maximal monotone graph as in the previous section; moreover, let α be a positive real number. We deal with the following problem.

Problem 5.1. To find a couple (u, s) such that u ∈ K and s ∈ H^{1}(Q) ∩ L^{∞}(Q), s(x, 0) = s^{0}(x) a.e. in Ω, and

∫∫_{Q} [s_{t}(u − v) + k(s)∇(u + ρgz) · ∇(u − v)]dxdt ≤ 0 ∀ v ∈ K  (5.1)
αs_{t} + G(s, s^{0}) + β(s) ∋ u a.e. in Q.  (5.2)

The variational inequality (5.1) is a weak formulation of equation (2.2) coupled with the Signorini condition (2.7). Inclusion (5.2) is the relaxed saturation versus pressure constitutive relation (see Section 3). In the sequel, for simplicity of notation, we shall omit the constant factor ρg in front of z.

Our aim is to prove an existence result for Problem 5.1.

We assume that

G(·, s^{0}) : C^{0}([0, T]) → C^{0}([0, T])

is bounded and Lipschitz continuous and that G+β is the inverse of a hysteresis operator F which satisfies the property (omitting the initial state among the entries of F)

∫_{0}^{T} [v(·, t) − v(·, 0)]\frac{∂F(v)}{∂t} dt ≥ 0  a.e. in Ω, ∀ v : Q → ℝ suitably regular. (5.3)

In some physical settings, an inequality like (5.3) accounts for energy dissipation along counterclockwise hysteresis cycles. See Visintin [22: p. 271] for regularity specification. However, in the sequel we shall suppose that the following discrete variant of (5.3) holds.

Let v ∈ L^{2}(Ω; C^{0}([0, T])) and {t_{n}}_{0≤n≤m} be a partition of [0, T] of uniform step-size h. Let v_{m} denote the linear time interpolate of the nodal values v(·, t_{n}) and, omitting the initial state among the entries of F, let us define w_{m}^{n} = [F(v_{m})](t_{n}) a.e. in Ω for every n. Then

h ∑_{n=1}^{m} \frac{w_{m}^{n} − w_{m}^{n−1}}{h} v(x, t_{n}) ≥ r  a.e. in Ω  (5.4)

where r is a real constant independent of the partition. In Section 6, we shall see that the Preisach hysteresis operator satisfies (5.3) - (5.4); that operator seems also appropriate to represent hysteresis in filtration phenomena.

Finally, given ŝ > 0 as in the previous section, we suppose that there exists ŝ < s̃ ≤ 1 such that

s̃ ≤ s^{0} ≤ 1  a.e. in Ω  (5.5)

and that there exists Γ̃ ⊂ Γ with positive bidimensional measure such that

P (x, t) > 0  a.e. on Γ̃, ∀ t ∈ [0, T].  (5.6)

Theorem 5.1. Let us suppose that the hysteresis operator G is bounded and Lipschitz continuous in C^{0}([0, T]), that G + β is the inverse of a hysteresis operator F satisfying (5.4), and finally that (5.5) - (5.6) hold. Then Problem 5.1 admits a solution.

We shall first prove Theorem 5.1 assuming that the quantity k(s) is larger than a positive constant. This will allow us to make use of a “maximum principle”. Then, in Lemma 5.2, we shall prove that the assumption k(s) ≥ c > 0 is correct.

Proof of Theorem 5.1. Approximation. Let us fix any m ∈ ℕ. We shall use an implicit time discretization of step-size h = \frac{T}{m}. Let P_{m}^{n} be the piecewise constant approximation of P, and ψ_{m}^{n} be the discretization of the hysteresis operator G, which we defined in Section 4. At every step n we suppose that s_{m}^{n−1} is known, and consider the problem of finding (u_{m}^{n}, s_{m}^{n}) such that s_{m}^{n} ∈ L^{∞}(Ω) and

\begin{matrix}u_{m}^{n} ∈ K_{m}^{n} := \{v ∈ H^{1}(Ω) : (γ_{0}v)^{+} = P_{m}^{n} on Γ\} and ∀ v ∈ K_{m}^{n} \\ ∫_{Ω} [\frac{s_{m}^{n} − s_{m}^{n−1}}{h} (u_{m}^{n} − v) + k(s_{m}^{n−1})∇(u_{m}^{n} + z) · ∇(u_{m}^{n} − v)]dx ≤ 0 \\ α\frac{s_{m}^{n} − s_{m}^{n−1}}{h} + ψ_{m}^{n}(s_{m}^{n}) + β(s_{m}^{n}) ∋ u_{m}^{n}. \end{matrix} \}  (5.7)

We solve problem (5.7) by a fixed-point procedure. In (5.7)_{3} we replace u_{m}^{n} by a function f ∈ L^{2}(Ω), and note that for every such f there exists a unique solution s_{f} of the equation (see Section 4). Moreover, since k(s_{m}^{n−1}) ≥ c > 0 (see Lemma 5.2 below), for every such s_{f} replacing s_{m}^{n}, there exists a unique solution u_{f} of (5.7)_{2}. The latter coincides with the unique minimizer in K_{m}^{n} of the following continuous strictly convex and coercive functional J_{f} which is defined on H^{1}(Ω), equipped with the equivalent norm ∥u∥^{2} = ∫_{Ω} |∇u|^{2}dx + ∫Γ̃ |γ_{0}u|^{2}dσ:

J_{f}(v) = ∫_{Ω} [\frac{k(s_{m}^{n−1})}{2}|∇(v + z)|^{2} + \frac{s_{f} − s_{m}^{n−1}}{h} v]dx.  (5.8)

The coercivity of J_{f} on K_{m}^{n} comes from the Poincaré inequality in view of (5.6) and definition (5.7)_{1}. Thus we have a function A : f ↦ u_{f}. To prove that (5.7) has a solution, it suffices to show that A has a fixed point in K_{m}^{n}. First of all, since u_{f} is a minimizer of J_{f}, it is easy to see that A(K_{m}^{n}) is bounded in H^{1}(Ω). Indeed, since k is bounded and s_{f} ∈ [0, 1] a.e. in Ω for every f ∈ L^{2}(Ω), we have that J_{f}(u_{f}) ≤ C with C independent on f; in particular, when f ∈ K_{m}^{n}, by coercivity we obtain uniform boundedness of u_{f} in H^{1}(Ω). Hence A(K_{m}^{n}) has compact closure in L^{2}(Ω).

Now we have to prove the continuity of A in L^{2}(Ω). Let f_{ℓ} converge to f in L^{2}(Ω). By the Lipschitz property of the operator F (cf. Proposition 4.2), we get that s_{ℓ} = s_{f_{ℓ}} converges to s_{f} in L^{2}(Ω). Moreover, the sequence u_{ℓ} = u_{f_{ℓ}} is bounded in H^{1}(Ω); hence (possibly extracting a subsequence) it converges to a function u weakly in H^{1}(Ω) and strongly in L^{2}(Ω). For every ℓ we have J_{ℓ}(u_{ℓ}) := J_{f_{ℓ}}(u_{ℓ}) ≤ J_{ℓ}(u_{f}) and then, using also the semicontinuity of ∫_{Ω} |∇(u_{ℓ} + z)|^{2}dx, we get

J_{f}(u) ≤ lim inf_{ℓ→+∞} J_{ℓ}(u_{ℓ}) ≤ J_{f}(u_{f}).

This yields u = u_{f}, and thus the continuity is proved. By Schauder’s theorem, there exists u_{m}^{n} ∈ K_{m}^{n} such that A(u_{m}^{n}) = u_{m}^{n}. Finally, let us set s_{m}^{n} = s_{u_{m}^{n}}.

A priori estimates. We use the same notation as in the previous section to denote linear time interpolations and piecewise constant functions. Let us assume that h < 1 and take in (5.7)_{1} v = (1 − h)u_{m}^{n} + hP_{m}^{n} ∈ K_{m}^{n}. We obtain

∫_{Ω} [(s_{m}^{n} − s_{m}^{n−1})(u_{m}^{n} − P_{m}^{n}) + hk(s_{m}^{n−1})∇(u_{m}^{n} + z) · ∇(u_{m}^{n} − P_{m}^{n})] ≤ 0.  (5.9)

Let us sum from 1 to n. Note that

\begin{matrix}∑_{ℓ=1}^{n} ∫_{Ω}(s_{m}^{ℓ} − s_{m}^{ℓ−1})u_{m}^{ℓ}dx \\ = h ∑_{ℓ=1}^{n} ∫_{Ω} [\frac{s_{m}^{ℓ} − s_{m}^{ℓ−1}}{h} (u_{m}^{ℓ} − α\frac{s_{m}^{ℓ} − s_{m}^{ℓ−1}}{h} ) + α\frac{(s_{m}^{ℓ} − s_{m}^{ℓ−1})^{2}}{h^{2}} ]dx. \end{matrix}  (5.10)

By the definition of the discretized hysteresis operator and by the discrete property (5.4), we have

h ∑_{ℓ=1}^{n} ∫_{Ω} \frac{s_{m}^{ℓ} − s_{m}^{ℓ−1}}{h} (u_{m}^{ℓ} − α\frac{s_{m}^{ℓ} − s_{m}^{ℓ−1}}{h} )dx ≥ r

where r is a real number independent from n and h.

Moreover, by discrete partial integration, for a suitable constant C independent on n and h, we get

\begin{matrix}|∑_{ℓ=1}^{n} ∫_{Ω}(s_{m}^{ℓ} − s_{m}^{ℓ−1})P_{m}^{ℓ}dx| \\ ≤ h ∑_{ℓ=1}^{n−1} ∫_{Ω} |s_{m}^{ℓ} \frac{P_{m}^{ℓ+1} − P_{m}^{ℓ}}{h} |dx + h ∫_{Ω} |s_{m}^{n}P_{m}^{n} − s^{0}P_{m}^{1}| dx \\ ≤ C. \end{matrix}  (5.11)

This quantity is uniformly bounded because P_{m} (the piecewise linear in time approximate of P) converges to P in H^{1}(0, T; L^{2}(Ω)). Using (5.11) and the fact that k(s_{m}^{n−1}) ≥ c > 0 for every n (see Lemma 5.2 below), via standard techniques we get

h ∑_{ℓ=1}^{n} ∫_{Ω} |∇u_{m}^{n}|^{2}dx ≤ C  (5.12)

where the constant C is independent on n and h. Hence (denoting the piecewise constant in time and the piecewise linear in time approximates as in the previous section)

∥\overline{u}_{m}∥_{L^{2}(0,T;H^{1}(Ω))} ≤ C,  independent on m.  (5.13)

Moreover, integrating inequality (4.6) in space and in time and using (5.13) we obtain ∥s_{m}∥_{H^{1}(0,T;L^{2}(Ω))} ≤ C, independent on h. Finally (see the end of Section 4), s_{m}^{n} ∈ H^{1}(Ω) for every m and n. Hence, s_{m} ∈ L^{2}(0, T; H^{1}(Ω)) and ∥s_{m}∥_{L^{2}(0,T;H^{1}(Ω))} ≤ C for every m. Therefore

∥s_{m}∥_{H^{1}(Q)} ≤ C  independent on m.  (5.14)

Passage to the limit. By (5.13) - (5.14) there exist two functions u ∈ L^{2}(0, T; H^{1}(Ω)) and s ∈ H^{1}(0, T; L^{2}(Ω)) such that, possibly extracting subsequences,

\overline{u}_{m} → u weakly in L^{2}(0, T; H^{1}(Ω))
s_{m} → s weakly in H^{1}(Q), strongly in C^{0}([0, T]; L^{2}(Ω));

(5.15)

the latter convergence follows from the compactness of the injection of H^{1}(Q) in the space C^{0}([0, T]; L^{2}(Ω)). By a procedure similar to (4.8) - (4.9), we get that the couple (u, s) solves (5.2).

Let us now take v as in (5.1), in particular assume that v ∈ C^{0}([0, T]; H^{1}(Ω)); this is not restrictive by a density argument. By (5.7) we get

∫∫_{Q} [\frac{∂s_{m}}{∂t} (\overline{u}_{m} − \overline{v}_{m}) + k(\underline{s}_{m})∇(\overline{u}_{m} + z) · ∇(\overline{u}_{m} − \overline{v}_{m})]dxdt ≤ 0  (5.16)

where \underline{s}_{m}(x, t) := \overline{s}_{m}(x, t − h) a.e. in Q (with the position \overline{s}_{m}(x, δ) := s^{0}(x) for all δ ∈ [−h, 0] a.e. in Ω). We claim that, passing to the inferior limit as m → +∞ in (5.16), we get (5.1) (with 1 in place of ρg). By (5.15), possibly extracting a further subsequence, we can suppose that \underline{s}_{m} → s a.e. in Q, hence k(\underline{s}_{m}) → k(s) a.e. in Q. Hence

\begin{matrix}lim inf_{m→+∞} ∫∫_{Q} k(\underline{s}_{m})∇(\overline{u}_{m} + z) · ∇(\overline{u}_{m} − \overline{v}_{m})dxdt \\ ≥ ∫∫_{Q} k(s)∇(u + z) · ∇(u − v)dxdt. \end{matrix}  (5.17)

By (5.7)_{3}, considering piecewise constant interpolation, we have

\overline{u}_{m} = α\frac{∂s_{m}}{∂t} + \overline{ψ}_{m} + ξ_{m}  (5.18)

where ξ_{m} ∈ β(\overline{s}_{m}). Hence, inserting (5.18) into the parabolic term of (5.16), we claim that we can pass to the inferior limit. In fact, \overline{ψ}_{m} → G(s, s^{0}) strongly in L^{2}(Q) and

lim inf_{m→+∞} ∫∫_{Q} (\frac{∂s_{m}}{∂t} )^{2}dxdt ≥ ∫∫_{Q} (\frac{∂s}{∂t})^{2}dxdt.  (5.19)

Moreover, for any finite partition {t_{i}} of [0, T] and for any selection ξ ∈ β(s), ξ ∈ L^{2}(Q), for every i there holds

\begin{aligned}∫_{Ω}[s(x, t_{i}) − s(x, t_{i−1})]ξ(x, t_{i−1}) dx & ≤ ∫_{Ω}[b(s(x, t_{i})) − b(s(x, t_{i−1}))]dx \\ & ≤ ∫_{Ω}[s(x, t_{i}) − s(x, t_{i−1})]ξ(x, t_{i}) dx. \end{aligned}  (5.20)

Passing to the limit in (5.20) as max{t_{i} − t_{i−1}} → 0, we get

∫∫_{Q} \frac{∂s}{∂t}ξdxdt = ∫_{Ω}[b(s(x, T)) − b(s^{0}(x))] dx.  (5.21)

Using (5.21) we finally obtain the inequality

lim inf_{m→+∞} ∫∫_{Q} \frac{∂s_{m}}{∂t} ξ_{m}dxdt ≥ lim inf_{m→+∞} ∑_{n=1}^{m} ∫_{Ω} [b(s_{m}^{n}(x)) − b(s_{m}^{n−1}(x))]dx

= lim inf_{m→+∞} ∫_{Ω} [b(s_{m}(x, T)) − b(s^{0}(x))]dx

≥ ∫_{Ω} [b(s(x, T)) − b(s^{0}(x))]dx

= ∫∫_{Q} \frac{∂s}{∂t}ξ dxdt

from which the proof is complete

Next, we state and prove a lemma which we already used to prove that k(s) is larger than a positive constant.

Lemma 5.2. Under the assumptions of Theorem 5.1, there exists ŝ < s_{min} < 1 independent of n (and m) such that s_{m}(x) ≥ s_{min} a.e. in Ω.

Proof. Let us define Z = sup_{Ω} z, where the function z is the vertical coordinate (we still omit the constant ρg in front of z). Let F̃ = (G + β)^{−1} represent the hysteretic relation between u and s (see Figure 4). By hypothesis, for a suitable C > 0, F̃ does not present hysteresis in (−∞, −C] and we may assume F̃(−C) < s̃.

By (5.7)_{3}, for a.e. x ∈ Ω, if (u_{m}^{n}(x), s_{m}^{n}(x)) stays above the graph F̃, then the term s_{m}^{n}(x) − s_{m}^{n−1}(x) is negative (the time derivative of s forces the couple to go towards the graph); on the contrary, if (u_{m}^{n}(x), s_{m}^{n}(x)) stays under the graph, then s_{m}^{n}(x) − s_{m}^{n−1}(x) is positive.

Let us consider the first step n = 1, and note that k(s^{0}) ≥ c > 0 a.e. in Ω. Then we take in (5.7)_{2} v = u_{m}^{1} + (u_{m}^{1} + z + C)^{−} ∈ K_{m}^{n} and obtain

∫_{Ω} [ − \frac{s_{m}^{1} − s^{0}}{h} (u_{m}^{1} + z + C)^{−} + k(s^{0})|∇(u_{m}^{1} + z + C)^{−}|^{2}]dx ≤ 0.

Again, for a.e. x ∈ Ω, if (u_{m}^{1}(x), s_{m}^{1}(x)) stays above the graph F̃, then the first term in the integral is nonnegative; on the other hand, if (u_{m}^{1}(x), s_{m}^{1}(x)) stays under the graph F̃, then the first term in the integral is still non-negative since it does not vanish only for u_{m}^{1}(x) < −z − C < −C and in that case we necessarily get s_{m}^{1}(x) < F̃(−C) ≤ s^{0}. Hence, we obtain

∇[(u_{m}^{1} + z + C)^{−}] = 0 a.e. in Ω.

Since (u_{m}^{1} + z + C)^{−} = 0 on a fixed part of ∂Ω with positive measure (see (5.6) and the definition of K), we get u_{m}^{1} + z ≥ −C a.e. in Ω and in particular u_{m}^{1} ≥ −C − Z a.e. in Ω. Now, we define s_{min} = F̃(−C − Z) and claim that s_{m}^{1} ≥ s_{min} a.e. in Ω. In fact, by (5.7)_{3}, for a.e. x ∈ Ω, s_{m}^{1}(x) < s_{min} implies that (u_{m}^{1}(x), s_{m}^{1}(x)) stays under the graph F̃ and hence s^{0}(x) < s_{m}^{1}(x) which is a contradiction.

Let us now suppose that u_{m}^{ℓ} ≥ −C − z and s_{m}^{ℓ} ≥ s_{min} a.e. in Ω for every 1 ≤ ℓ ≤ n−1. By induction, we show that the same inequalities hold for u_{m}^{n} and s_{m}^{n}, respectively. Let us take x_{0} ∈ Ω such that all the functions are defined for all 1 ≤ ℓ ≤ n and satisfy the inductive hypotheses and that (5.7)_{3} holds for all 1 ≤ ℓ ≤ n. By contradiction, let us suppose that u_{m}^{n}(x_{0}) < −C − z(x_{0}) and that (u_{m}^{n}(x_{0}), s_{m}^{n}(x_{0})) stays under the graph. In particular, it follows that s_{m}^{n−1}(x_{0}) < s_{m}^{n}(x_{0}) < s̃. By inductive hypothesis, u_{m}^{n−1}(x_{0}) ≥ −C − z(x_{0}) > u_{m}^{n}(x_{0}) and then also (u_{m}^{n−1}(x_{0}), s_{m}^{n−1}(x_{0})) stays under the graph. Hence, s_{m}^{n−2}(x_{0}) < s̃. Acting in this way, finally we obtain that s^{0}(x_{0}) < s̃ which is a contradiction. Then we can conclude that (u_{m}^{n}(x), s_{m}^{n}(x)) stays above the graph at every almost points x ∈ Ω such that u_{m}^{n}(x) < −C − z(x). As in the first step, taking the test function v = u_{m}^{n} + (u_{m}^{n} + z + C)^{−} in (5.7), we obtain u_{m}^{n} ≥ −C − z a.e. in Ω and in particular u_{m}^{n} ≥ −C − Z a.e. in Ω. Arguing as before, we can then conclude that, at almost every point, s_{m}^{n} ≥ s_{min} otherwise we should obtain s^{0} < s_{min} < s̃ which is again a contradiction

Remark 5.3. Some engineering papers (see, for instance, Poulovassilis and Tzimas [20], Kacimov and Yakimov [10], and the references therein) account for hysteresis even in the relation between hydraulic conductivity and saturation, which in this paper is represented as in Figure 3 (without hysteresis). These hysteretic effects are however much less evident than the ones in the s versus u relation.

Our result works in this case too. Indeed, let us call K(·, s^{0}) the hysteresis operator representing the relation between k and s, and suppose that it is bounded and Lipschitz continuous in C^{0}([0, T]) (again we consider x just as a parameter). We replace in (5.8) k(s_{m}^{n−1}) by the approximation k_{m}^{n−1}(s_{m}^{n−1}) of the operator K (constructed as the approximation of G in Section 4). In particular, note that k_{m}^{n−1}(s_{m}^{n−1}) is to be considered as known in (5.8). Hence we again obtain (5.15) and (5.17) (with k(\underline{s}_{m}) and k(s) replaced by k_{m}^{n−1}(s_{m}^{n−1}) and K(s, s^{0}), respectively), from which the conclusion follows.

Remark 5.4. A more general model should consider the hydraulic conductivity as a tensor depending on the saturation s and on the point x (the model we studied here corresponds to the case of isotropic material and independence on x). However, our result can be easily extended to the case where the hydraulic conductivity is given by k̃(s, x) = a(x)k(s) with k as in Figure 3 and a(x) uniformly (on x) strictly positive definite matrix.

6. The hysteretic component of the constitutive relation

As we pointed out, the engineering literature seems to support the use of the Preisach model in the description of the saturation versus pressure constitutive relation (in particular we refer to what they call domain theory of hysteresis). This is not surprising as this model has been applied to several phenomena, after it was proposed in the 1930s by the physicist F. Preisach to represent scalar ferromagnetic hysteresis.

The idea at the basis of this model is simple and appealing: a hysteresis loop is seen as the superposition of a family of rectangular loops, just as a real function can be represented as an average of shifted and weighted jumps: for any f ∈ W^{1,1}(ℝ),

f(x) = lim_{s→−∞} f(s) + ∫_{ℝ} H(x − s)f′(s)ds

(here H is the Heaviside function: H(s) = 0 if s ≤ 0 and H(s) = 1 if s > 0).

Delayed Relays. Let us denote by C_{r}^{0}([0, T]) the space of functions [0, T] → ℝ which are continuous on the right in [0, T [. For any pair ρ = (ρ_{1}, ρ_{2}) ∈ ℝ^{2} (ρ_{1} < ρ_{2}), we introduce the delayed relay operator

h_{ρ} : C^{0}([0, T]) × {−1, 1} → BV (0, T) ∩ C_{r}^{0}([0, T]).

For any u ∈ C^{0}([0, T]) and any ξ, ξ = −1 or ξ = 1, the function z = h_{ρ}(u, ξ) : [0, T] → {−1, 1} is defined by

z(0) = \{ \begin{matrix}−1 & if u(0) ≤ ρ_{1} \\ ξ & if ρ_{1} < u(0) < ρ_{2} \\ 1 & if u(0) ≥ ρ_{2} \end{matrix}  (6.1)

and for any t ∈ (0, T], setting X_{t} = {τ ∈ (0, t] : u(τ) = ρ_{1} or u(τ) = ρ_{2}},

z(t) = \{ \begin{matrix}z(0) & if X_{t} = ∅ \\ −1 & if X_{t} ≠ ∅ and u(max X_{t}) = ρ_{1} \\ 1 & if X_{t} ≠ ∅ and u(max X_{t}) = ρ_{2}. \end{matrix}  (6.2)

Thus z is uniquely defined in [0, T]. For instance, let u(0) < ρ_{1}. Then z(0) = −1, and z(t) = −1 as long as u(t) < ρ_{2}; if at some instant u reaches ρ_{2}, then z jumps up to 1, where it remains as long as u(t) > ρ_{1}; if later u reaches ρ_{1}, then z jumps down to −1; and so on (cf. Figure 6).

Figure 6: Delayed relay

The Preisach Model. The thresholds of delayed relay operators form the so called Preisach (half) plane

P = {ρ = (ρ_{1}, ρ_{2}) ∈ ℝ^{2} : ρ_{1} < ρ_{2}} .  (6.3)

We denote by R the family of Borel measurable functions P → {−1, 1}, and by ξ = {ξ_{ρ}} a generic element of R, which we intend to represent the initial configurations of all the delayed relays. We fix a finite (signed) Borel measure μ over P, and introduce the corresponding Preisach operator

\begin{matrix}H_{μ} : C^{0}([0, T]) × R → L^{∞}(0, T) ∩ C_{r}^{0}([0, T]) \\ [H_{μ}(u, ξ)](t) = ∫_{P}[h_{ρ}(u, ξ_{ρ})](t) dμ(ρ) ∀ t ∈ [0, T]. \end{matrix} \}  (6.4)

If μ ≥ 0, it is not difficult to see that all the hysteresis branches (non only the exterior loop, but also the interior loops) are non-decreasing (so-called piecewise monotonicity). In the following statement, |μ| denotes the total variation of the measure μ.

Theorem 6.1 (Continuity). Let μ be a finite Borel measure over P and ξ ∈ R. Then H_{μ}(u, ξ) ∈ C^{0}([0, T]) for any u ∈ C^{0}([0, T]) if and only if |μ|([r, +∞)) × {r}) = |μ|({r} × [r, +∞)) = 0 for any r ∈ ℝ. Moreover, if the latter condition holds, then H_{μ}(·, ξ) is strongly continuous in C^{0}([0, T]).

For the proof we refer to Visintin [22: p. 113].

Under appropriate conditions on the measure μ, the operator H_{μ}(·, ξ) is continuous, or uniformly continuous, or Lipschitz continuous in C^{0}([0, T]), or operates in the Sobolev spaces W^{1,p}(0, T) (1 ≤ p ≤ +∞), or in the Hölder spaces C^{0,ν}([0, T]) (0 < ν ≤ 1), or in C^{0}([0, T]) ∩ BV (0, T) (see Brokate and Visintin [7]).

The two following results allow us to apply the results of this paper to the case in which the hysteresis relation is represented by the Preisach model. The first one was proved by Brokate and Visintin [7: Theorem 5.14].

Theorem 6.2 (Lipschitz Continuity of the Inverse). Let μ be a finite non-negative Borel measure over P, ξ ∈ R and a > 0. Then aI + H_{μ}(·, ξ) : C^{0}([0, T]) → C^{0}([0, T]) is invertible, and its inverse operator is Lipschitz continuous in C^{0}([0, T]).

The following statement is easily checked, as it holds for each delayed relay.

Proposition 6.3. Let μ be a finite non-negative Borel measure over P and ξ ∈ R. Then H_{μ}(·, ξ) : C^{0}([0, T]) → C^{0}([0, T]) fulfills inequalities (5.3) and (5.4).

Remark 6.4. The case of relays switching between two values different from +1 and −1 corresponds to an affine transformation of the measure μ. Hence (see Section 3) for a suitable non-negative measure μ vanishing out of the triangle of the Preisach plane Δ = {u_{1} ≤ ρ_{1} ≤ ρ_{2} ≤ u_{2}}, and for suitable nonnegative values h_{1} and h_{2} for the relays such that s_{1} = h_{1}μ(Δ) and s_{2} = h_{2}μ(Δ), we can invert the operator aI +H_{μ}(·, ξ) between the spaces C^{0}([0, T]; [u_{1}, u_{2}]) and C^{0}([0, T]; [s_{1}, s_{2}]. According to Theorem 6.2, the inverse operator is Lipschitz continuous. This allowed us in Sections 3 and 4 to suppose that the inverse hysteresis operator G̃ = F̃^{−1} is of the form G̃ = G + β, where G is a Lipschitz continuous hysteresis operator and β is a maximal monotone graph. Indeed, this holds whenever the hysteresis branches of F̃ are confined to a bounded subset of ℝ^{2} and have a minimum slope a > 0. In this case (see (3.6)) we can write F̃ = F ∘ τ + ϕ where ϕ is a maximal monotone graph, τ is a truncation on ℝ and F = aI + H, with H a hysteresis operator between the spaces C^{0}([0, T]; [u_{1}, u_{2}]) and C^{0}([0, T]; [s_{1}, s_{2}]).

Remark 6.5. In the previous sections we did not display the initial internal variable ξ of the Preisach operator, but we only displayed the initial output’s state s^{0}. This should be regarded as a “contraction” of a more precise formulation which should contain also the initial state of each relay, represented by the function ξ. By displaying only the initial variable s^{0}, we simplified the notation without restricting the generality of our results. More precisely, if we suppose that the initial internal variable ξ is a datum of the problem, then we can regard the hysteresis operator as depending only on the initial output’s state (besides the input function u). In particular (see (3.5) and Remark 6.4), when the initial state of the output is equal to s_{1} or s_{2}, all the relays are switched to h_{1} or h_{2}, respectively.

Hysteresis in Filtration through Porous Media 997

References

[1] Alt, H. W.: Strömungen durch inhomogene poröse Medien mit freiem Rand. J. Reine Angew. Math. 305 (1979), 89 – 115.

[2] Alt, H. W., Luckhaus, S. and A. Visintin: On nonstationary flow through porous media. Ann. Mat. Pura Appl. 136 (1984), 303 – 316.

[3] Bagagiolo, F.: Feedback control of a filtration problem with hysteresis. Adv. Math. Sci. Appl. 9 (1999), 787 – 804.

[4] Baiocchi, C.: Su un problema di frontiera libera connesso a questioni di idraulica. Ann. Mat. Pura Appl. 92 (1972), 107 – 127.

[5] Bear, J.: Dynamic of Fluids in Porous Media. New York: American Elsevier 1972.

[6] Brokate, M. and J. Sprekels: Hysteresis and Phase Transitions. Berlin: Springer-Verlag 1996.

[7] Brokate, M. and A. Visintin: Properties of the Preisach model for hysteresis. J. Reine Angew. Math. 402 (1989), 1 – 40.

[8] Fredlund, D. G. and H. Rahardjo: Soil Mechanics for Unsaturated Soils. New York: Wiley 1993.

[9] Gilardi, G.: A new approach to evolution free boundary problems. Comm. Part. Diff. Equ. 4 (1979), 1099 – 1122.

[10] Kacimov, A. B. and N. D. Yakimov: Nonmonotonic moisture profile as a solution of Richards’ equation for soils with conductivity hysteresis. Adv. Water Resour. 21 (1998), 691 – 696.

[11] Koll, J. B. and J. C. Parker: Development and evaluation of closed-form expressions for hysteretic soil hydraulic properties. Water Resour. Res. 23 (1987), 105 – 114.

[12] Krasnosel’skiĭ, M. A. and A. V. Pokrovskiĭ: Systems with Hysteresis. Berlin: Springer-Verlag 1989 (Russian ed.: Moscow: Nauka 1983).

[13] Krejčí, P.: Convexity, Hysteresis and Dissipation in Hyperbolic Equations. Tokyo: Gakkotosho 1996.

[14] Mayergoyz, I. D.: Mathematical Models of Hysteresis. New York: Springer-Verlag 1991.

[15] Mualem, Y.: Theory of universal hysteretical properties of unsaturated porous media. In: Proc. Fort Collins Fluid Intern. Hydrol. Symp. 1977, Water Resour. Res. Public. 1979.

[16] Mualem, Y.: A modified dependent domain theory of hysteresis. Soil Sci. 137 (1984), 283 – 291.

[17] Otto, F.: L^{1}-contraction and uniqueness for unstationary saturated-unsaturated porous media flow. Adv. Math. Sci. Appl. 7 (1997), 537 – 553.

[18] Poulovassilis, A. and E. C. Childs: The hysteresis of pore water: the non-independence of domains. Soil Sci. 112 (1971), 301 – 312.

[19] Poulovassilis, A. and W. M. El-Ghamry: The dependent domain theory applied to scanning curves of any order in hysteretic soil water relationships. Soil Sci. 126 (1978), 1 – 8.

[20] Poulovassilis, A. and E. Tzimas: The hysteresis in the relationship between hydraulic conductivity and soil water content. Soil Sci. 120 (1975), 327 – 331.

[21] Torelli, A.: Su un problema a frontiera libera di evoluzione. Boll. Un. Mat. Ital. 11 (1975), 559 – 570.

[22] Visintin, A.: Differential Models of Hysteresis. Heidelberg: Springer-Verlag 1994.

998 F. Bagagiolo and A. Visintin Received 27.12.1999