Deformable porous media with degenerate hysteresis in gravity field

1.   Introduction

We pursue here the study started in ^[12,13] of degenerate diffusion in unsaturated porous media filled with liquid and gas. Here we additionally take into account the deformations of the solid matrix produced by the penetrating humidity, and the effects of gravity on fluid diffusion as a generalization of the Richards equation (see ^[24]). To be precise, our main modeling assumptions are the following:

${\rm (a)}$ The pressure-saturation relation exhibits hysteresis;

${\rm (b)}$ The solid skeleton is viscoelastic;

${\rm (c)}$ The permeability of the material depends on the moisture content;

${\rm (d)}$ A gravity-driven moisture flux takes part in the process.

Our model is inspired by the experimental evidence about hysteresis in hydrogeology in, e.g., ^[2,16,20]. Figure 1 is taken from ^[20] and shows rate-independent hysteresis between the logarithm soil suction $\psi$, which is a decreasing function of the pressure, and the volumetric water content $\theta$ which, up to a linear transformation, can be identified with the saturation. The main characteristic of the hysteresis relation in Figure 1 is that at turning points where the pressure changes the orientation from decreasing to increasing or vice versa, the starting slope is horizontal, which makes the diffusion problem degenerate.

Mathematical investigation of diffusion problems with hysteresis goes back to A. Visintin's pioneering monograph ^[25] presenting basic methods of solving PDEs with hysteresis. However, two problems have remained unsolved until recently: the degeneracy at turning points and the dependence of the diffusion coefficient on the saturation. We have proved in ^[12] that the degeneracy in the case of constant permeability can be overcome using a method developed there to convexify the hysteresis relation. In ^[13], the convexification combined with a compact embedding theorem for anisotropic Orlicz spaces was shown to give a sufficient argument for the solvability also in the case of saturation-dependent permeability. Let us also mention some previous results on hysteresis dependence of the permeability coefficient in the non-degenerate case under additional regularization in time or space ^[4,5,9,26]. To our knowledge, only E. El Behi-Gornostaeva addressed in her thesis ^[8] the full problem of hysteresis-dependent saturation without regularization in the non-degenerate case and proposed a method for the existence proof.

In this paper we extend the problem to the case of deformable porous media under gravity effects and full degeneracy of the saturation dependence. The process is driven by the mass conservation principle and the quasistatic mechanical equilibrium equation in the small deformations regime. We further simplify the system in order to make it mathematically tractable by assuming that shear stresses are negligible. The gravity term and the viscoelasticity of the deformable solid skeleton represent a new degree of complexity related to the fact that no obvious a priori upper bound for the solutions is available here. We show that such an estimate can be obtained using a hysteresis variant of the Moser iteration technique under an additional assumption on the admissible Preisach operators $G$, namely, that the Preisach density $\phi(x, r, v)$ in (3.6) of $G$ decays sufficiently slowly at infinity, see Hypothesis 3.6. The proof of the existence theorem is then carried out by a convexity and compactness argument developed in ^[13].

The structure of the present paper is the following. In Section 2 we set up the mathematical model for the phenomenon and show that it is consistent with the physical principles of mass and energy conservation. In Section 3 we recall the definitions of the main concepts including convexifiable Preisach operators, and state the main existence Theorem 3.7. In Section 4 we propose a time discretization scheme with time step $\tau > 0$ and derive estimates independent of $\tau$. A uniform upper bound for the time-discrete approximation is established in Section 5 via Moser-type iterations. In Section 6 we show that, similarly to ^[13], we can derive the convexity estimate (6.26) for the time derivative of the pressure. As shown in the last Section 7, this is sufficient to let the time step $\tau$ tend to $0$ and to prove that the limit is the desired solution.

2.   Derivation of the model

The physical quantities which appear in the model are listed below in Table 1.

The fluid diffusion in a domain $\Omega \subset {\mathbb{R}}^3$ filled with a deformable unsaturated porous medium is driven by the liquid mass conservation principle

with constant mass density $\rho$, and by the mechanical equilibrium equation for the solid

where ${\boldsymbol{n}}(x)$ denotes the unit outward normal vector to $\Omega$ at a point $x \in \partial\Omega$. We write (2.2) in variational form

where $\nabla_s$ denotes the symmetric gradient.

We assume that the deformations are small, so that ${\mathrm{\, div\, }} {\boldsymbol{u}}$ is the relative local volume change, and that shear stresses can be neglected, that is,

where $\mu > 0$ is a constant bulk elasticity modulus, $\gamma > 0$ is a bulk viscosity modulus, and ${\boldsymbol{\delta}} = (\delta_{ij})$, $i, j = 1, 2, 3$ is the Kronecker symbol. Then (2.3) can be written as

Let $h \in L^2(\Omega)$ be arbitrary, and let $W$ be the solution to the problem

Putting ${\boldsymbol{\xi}} = \nabla W$ in (2.5), we get the mechanical equilibrium equation in the form

a.e. in $\Omega$.

Assume now that the relative liquid mass flux ${\boldsymbol{q}} - \rho {\boldsymbol{u}}_t$ is proportional to the liquid pressure gradient $\nabla p$, that is,

with proportionality factor $\kappa = \kappa(x, \theta) > 0$ depending on $x \in \Omega$ and $\theta$. The liquid pressure $p$ can be decomposed into two components

where $p_c$ is the capillary pressure and $p_h$ is the hydrostatic pressure. For the hydrostatic pressure we assume the classical relation

where ${\boldsymbol{\nu}}$ is the unit vector in the gravity direction and $x_0$ is a referential point.

We introduce the dimensionless normalized capillary pressure

and follow the modeling hypotheses of ^[2,24] which consists in representing the $p_c \mapsto \theta$ hysteresis relation by the formula

where $G$ is a Preisach operator defined below in Section 3.

On the boundary $\partial\Omega$ of $\Omega$ we assume that the unit outward normal ${\boldsymbol{n}}(x)$ is defined almost everywhere and that the normal component of the relative flux ${\boldsymbol{q}} - \rho {\boldsymbol{u}}_t$ in (2.8) is proportional to the difference of pressures $p$ inside and $p^*$ outside the body, that is,

The weak formulation of the mass balance Eq (2.1) then reads

for every test function ${y} \in W^{1, 2}(\Omega) \cap L^\infty(\Omega)$.

Assume for the moment that there is no mass exchange with the exterior corresponding to the choice $b^*(x) = 0$ in (2.13), that is,

on $\partial\Omega$. Putting ${y} = 1$ in (2.14), we get

The interpretation of (2.16) is mass conservation. The term ${\int_{\Omega}} {\mathrm{\, div\, }} {\boldsymbol{u}}(x, t){\, \mathrm{d}} x$ describes the evolution of the volume represented by $\Omega$, and $\rho{\int_{\Omega}} (\theta + {\mathrm{\, div\, }} {\boldsymbol{u}})(x, t){\, \mathrm{d}} x$ is the total water mass in $\Omega$ at time $t$. Naturally, if the volume increases and mass is conserved, then the saturation decreases and vice versa.

Consider now the energy balance. Under the boundary condition (2.15), no power is supplied to the system, and the total energy of the system should therefore decrease because all phenomena like viscosity, diffusion, and hysteresis dissipate energy. We put ${y} = p/\rho$ in (2.14) and obtain using (2.15) that

For a Preisach operator $G$ in (2.12) there exists a Preisach potential operator $V$ and a dissipation operator $D$ such that

Using (2.7) and (2.9)–(2.11), we can rewrite (2.17) in the form

The energetic interpretation of (2.19) is the following: The term

is the dissipation rate which is positive in agreement with the principles of thermodynamics, while

is the potential energy of the system which is therefore decreasing as expected.

Putting $v = {\mathrm{\, div\, }} {\boldsymbol{u}} - p_h/\mu$, from (2.7) and (2.9)–(2.14) we get the following system, consisting of a PDE coupled with an ODE, for the unknowns $u$ and $v$:

for every ${y} \in W^{1, 2}(\Omega) \cap L^\infty(\Omega)$, where we set $u^* = (p^*- p_h)/p_0$.

The existence proof for (2.22)-(2.23) is independent of the exact values of the physical constants. We therefore normalize all constants to $1$ for simplicity and consider the system

for every ${y} \in W^{1, 2}(\Omega)\cap L^\infty(\Omega)$, coupled with the initial conditions

3.   Statement of the problem

We study here Problem (2.24)–(2.26) in a bounded Lipschitzian domain $\Omega \subset {\mathbb{R}}^N$ and time interval $(0, T)$, where $N \in {\mathbb{N}}$ can be arbitrary. Recall first the definition of the Preisach operator introduced in ^[21] in the variational setting of ^[18].

Definition 3.1. Let $\lambda \in L^\infty(\Omega \times (0, \infty))$ be a given function with the following properties:

For a given $r > 0$, we call the play operator with threshold $r$ and initial memory $\lambda$ the mapping which with a given function $u \in L^1(\Omega; W^{1, 1}(0, T))$ associates the solution $\xi^r\in L^1(\Omega; W^{1, 1}(0, T))$ of the variational inequality

with initial condition

and we denote

Given a measurable function $\phi :\Omega\times(0, \infty)\times {\mathbb{R}} \to [0, \infty)$ and a constant $\bar G \in {\mathbb{R}}$, the Preisach operator $G$ is defined as a mapping $G: L^2(\Omega; W^{1, 1}(0, T))\to L^2(\Omega; W^{1, 1}(0, T))$ by the formula

The Preisach operator is said to be regular if the density function $\phi$ of $G$ in (3.6) belongs to $L^\infty(\Omega\times (0, \infty)\times {\mathbb{R}})$, and there exist constants $\phi_1, \bar{\phi} > 0$ and a decreasing function $\phi_0: {\mathbb{R}}\to {\mathbb{R}}$ such that for all $U > 0$, all $x, x_1, x_2 \in \Omega$, and a.e. $(r, v)\in (0, U) \times (-U, U)$ we have

In applications, the natural physical condition $\theta = G[u] \in [0, 1]$ is satisfied for each input function $u$ if and only if $\bar G \in (0, 1)$ and

for a.e. $x\in \Omega$. However, for the existence result in Theorem 3.7 only (3.7) and (3.8) are substantial.

Let us mention the following classical result (see [18, Proposition Ⅱ.3.11]).

Proposition 3.2. Let $G$ be a regular Preisach operator in the sense of Definition 3.1. Then it can be extended to a Lipschitz continuous mapping $G: L^p(\Omega; C[0, T]) \to L^p(\Omega; C[0, T])$ for every $p \in [1, \infty)$.

The Preisach operator is rate-independent. Hence, for input functions $u(x, t)$ which are monotone in a time interval $t\in (a, b)$, a regular Preisach operator $G$ can be represented by a superposition operator $G[u](x, t) = B(x, u(x, t))$ with an increasing function $u \mapsto B(x, u)$ called a Preisach branch. Indeed, the branches may be different at different points $x$ and different intervals $(a, b)$. The branches corresponding to increasing inputs are said to be ascending (the so-called wetting curves in the context of porous media), the branches corresponding to decreasing inputs are said to be descending (drying curves).

Definition 3.3. Let $U > 0$ be given. A Preisach operator is said to be uniformly counterclockwise convex on $[-U, U]$ if for all inputs $u$ such that $|u(x, t)|\le U$ a.e., all ascending branches are uniformly convex and all descending branches are uniformly concave in the sense that there exists $\beta > 0$ such that for every $u \in (-U, U)$ we have the implications

in the sense of distributions, see ^[12].

A regular Preisach operator $G$ is called convexifiable if for every $U > 0$ there exist a uniformly counterclockwise convex Preisach operator $P$ on $[-U, U]$, positive constants $g_*(U), g^*(U), \bar{g}(U)$, and a twice continuously differentiable mapping $g:[-U, U] \to [-U, U]$ such that

and $G = P\circ g$.

A typical example of a uniformly counterclockwise convex operator is the so-called Prandtl-Ishlinskii operator characterized by positive density functions $\phi(x, r)$ independent of $v$, see [18, Section 4.2]. Operators of the form $P\circ g$ with a Prandtl-Ishlinskii operator $P$ and an increasing function $g$ are often used in control engineering because of their explicit inversion formulas, see ^[3,6,19]. They are called the generalized Prandtl-Ishlinskii operators (GPI) and represent an important subclass of Preisach operators. Note also that for every Preisach operator $P$ and every Lipschitz continuous increasing function $g$, the superposition operator $G = P\circ g$ is also a Preisach operator, and there exists an explicit formula for its density, see [17, Proposition 2.3]. Another class of convexifiable Preisach operators is shown in [12, Proposition 1.3].

The technical hypotheses on the data in (2.24)-(2.25) can be stated as follows.

Hypothesis 3.4. The boundary permeability $b^*$ belongs to $L^\infty(\partial\Omega)$, and is such that $b^*(x)\ge 0$ a.e. and ${\int_{\partial\Omega}} b^*(x){\, \mathrm{d}} s(x) > 0$. The boundary source $u^*$ belongs to $L^\infty(\partial \Omega\times(0, T))$; in addition, $u^*_t \in L^2(\partial \Omega\times(0, T))$. The permeability $\kappa:\Omega \times {\mathbb{R}} \to {\mathbb{R}}$ is a bounded Lipschitz continuous function, more precisely, there exist constants $\bar\kappa > 0$, $\kappa_1 > \kappa_0 > 0$ such that for all $\theta, \theta_1, \theta_2\in {\mathbb{R}}$ and all $x, x_1, x_2 \in \Omega$ we have

Note that we can rewrite (2.24)-(2.25) for ${y} \in W^{1, 2}(\Omega)\cap L^\infty(\Omega)$ as

Therefore, taking into account (2.26), even a local solution to Problem (2.24)–(2.26) may fail to exist if for example $\lambda(x, r) \equiv 0$ and $\!{\mathrm{\, div\, }}\!\big(\kappa(x, G[u](x, 0))\left(\nabla u_0(x) + {\boldsymbol{\nu}}\right)\big)+v_0(x)-u_0(x) \ne 0$, and we need an initial memory compatibility condition which we state here following ^[12]. A more detailed discussion on this issue can be found in the introduction to ^[12].

Hypothesis 3.5. Let the initial memory $\lambda$ and the Preisach density function $\phi$ be as in Definition 3.1. The initial pressure $u_0$ belongs to $W^{2, \infty}(\Omega)$, the initial volume strain $v_0$ belongs to $L^\infty(\Omega)$, and there exist a constant $L > 0$ and a function $r_0 \in L^\infty(\Omega)$ such that, for $\Lambda > 0$ as in (3.1), ${\mathop{\mathrm{\, sup\, ess}}}_{x\in \Omega}|u_0(x)| \le \Lambda$ and for a.e. $x \in \Omega$ the following initial compatibility conditions hold:

Unlike ^[12], here we do not need to assume $\!{\mathrm{\, div\, }}\!\big(\kappa(x, \theta_0(x))\left(\nabla u_0(x) + {\boldsymbol{\nu}}\right)\big) \in L^\infty(\Omega)$ since it follows from the fact that $u_0 \in W^{2, \infty}(\Omega)$ together with assumptions (3.2), (3.8), and Hypothesis 3.4 on $\kappa$. Instead, in addition to the hypotheses of ^[12], we have to assume a polynomial decay of the Preisach density function at infinity that will allow us to perform a hysteresis variant of the Moser iterations in Section 5.

Hypothesis 3.6. The Preisach density $\phi :\Omega\times(0, \infty)\times {\mathbb{R}} \to [0, \infty)$ in (3.6) is such that

Condition (3.16) is indeed compatible with (3.9) if $m > 2$, as

so that (3.9) is satisfied if for example $\bar G = 1/2$ and $\phi_0 \le 1- 2/m$.

Our main existence result reads as follows.

Theorem 3.7. Let Hypotheses 3.4–3.6 hold, and let $G$ be a convexifiable Preisach operator in the sense of Definition 3.3. Then there exists a solution $(u, v)$ to Problem (2.24)–(2.26) such that $u, v \in L^\infty(\Omega\times (0, T))$, $\nabla u \in L^2(\Omega\times (0, T); {\mathbb{R}}^N)$, $u_t$ and $\theta_t = G[u]_t$ belong to the Orlicz space $L^{\Phi_{log}}(\Omega\times (0, T))$ generated by the function $\Phi_{log}(v) = v\log(1+v)$, and $v_t \in L^\infty(\Omega\times (0, T))$.

Basic properties of Orlicz spaces are summarized in [13, Section 5]. For a more comprehensive discussion, we refer the interested reader to the monographs ^[1,22].

4.   Time discretization

We proceed as in ^[12], choose a discretization parameter $n \in {\mathbb{N}}$, define the time step $\tau = T/n$, and replace system (2.24)-(2.25) with the following time-discrete counterpart for the unknowns $\{u_i, v_i: i = 1, \dots, n\}$:

for every test function ${y} \in W^{1, 2}(\Omega)$, where $u^*_i(x) = u^*(x, i\tau)$ for $i\in \{1, \dots, n\}$, $u_0$ and $v_0$ are as in (2.26). Here, the time-discrete Preisach operator $G[u]_i$ is defined for an input sequence $\{u_i : i\in{\mathbb{N}}\cup\{0\}\}$ by a formula of the form (3.6), namely,

where $\xi^r_i$ denotes the output of the time-discrete play operator

defined as the solution operator of the variational inequality

with a given initial condition

similarly as in (3.3)-(3.4). Note that the discrete variational inequality (4.5) can be interpreted as weak formulation of (3.3) for piecewise constant inputs in terms of the Kurzweil integral, and details can be found in [10, Section 2].

Proposition 4.1. System (4.1)-(4.2) with initial conditions $u_0$ and $v_0$ as in (2.26) admits at least one solution $\{u_i, v_i\} \subset W^{1, 2}(\Omega)$ for each $i\in \{1, \dots, n\}$.

Proof. We proceed by induction, and assume that the sequences $\{u_i, v_i\}$ have already been constructed for $i = 1, \dots, j-1$. Note that there is no hysteresis in the passage from $u_{j-1}$ to $u_j$, so that there exists a function $G_j: \Omega\times {\mathbb{R}} \to {\mathbb{R}}$, continuous and nondecreasing in the second variable, such that $G[u]_j = G_j(x, u_j)$. Furthermore, by (4.2) we have

Hence, (4.1)-(4.2) can be interpreted as a quasilinear elliptic equation for the unknown $u { :=} u_j$ of the form

for every $y \in W^{1, 2}(\Omega)$ with given functions $a_j \in L^2(\Omega)$.

Let $\{e_k: k\in {\mathbb{N}}\} \subset W^{1, 2}(\Omega)$ be an orthonormal basis of $L^2(\Omega)$. It is convenient to choose eigenfunctions of the following problem, which is compatible with boundary conditions (2.13), namely

with eigenvalues $0 < \mu_1 < \mu_2 \le \dots$. For each $m \in {\mathbb{N}}$ we look for coefficients $u_k$, $k = 1, \dots, m$, such that the function

satisfies the identity

for every $k = 1, \dots, m$. The existence of a solution to (4.9) can be proved in a classical way using the Brouwer degree theory. An introduction to topological methods for solving nonlinear partial differential equations can be found in [11, Chapter Ⅴ], but we will mainly refer to ^[14], where detailed proofs are also provided.

We begin by noting that we can define a continuous mapping $T: {\mathbb{R}}^m \to {\mathbb{R}}^m$ associating to $\boldsymbol{u} { :=} (u_1, \dots, u_m)$ the left-hand side of (4.9) for $k = 1, \dots, m$. Let $\gamma \in [0, 1]$, and consider the homotopy

Testing (4.10) by $u_k$ and summing up over $k = 1, \dots, m$, by Hypothesis 3.4 we obtain

Hence, by Young's inequality and again Hypothesis 3.4, we see that $T_\gamma(\boldsymbol{u})\cdot\boldsymbol{u} > 0$ outside the ball $B{^{(m)}}_R \subset {\mathbb{R}}^m$ with a sufficiently large radius $R$ independent of $\gamma$ and $m$. Therefore, the equation $T_\gamma(\boldsymbol{u}) = 0$ has no solution on $\partial B{^{(m)}}_R$, so that by [14, Theorem 3] the Brouwer degree $d(T_\gamma, B{^{(m)}}_R, 0)$ is a constant for $\gamma \in [0, 1]$. The mapping $T_0$ is linear, hence its degree is odd. In particular, $d(T_1, B{^{(m)}}_R, 0) = d(T_0, B{^{(m)}}_R, 0) \ne 0$, which implies (see [14, Theorem 2]) that the equation $T_1(\boldsymbol{u}) = 0$, equivalently, (4.9), has at least one solution in $B{^{(m)}}_R$.

Testing (4.9) by $u_k$, summing up over $k = 1, \dots, m$, and employing Hypothesis 3.4 we see that the sequence $\{u{^{(m)}}\}$ is uniformly bounded in $W^{1, 2}(\Omega)$, hence it is compact in $L^2(\Omega)$. Let $u$ be the limit of any convergent subsequence of $\{u{^{(m)}}\}$. Passing to the limit in (4.9) as $m\to \infty$, we check that (4.8) holds with $y = e_k$ for every $k \in {\mathbb{N}}$. Since the system $\{e_k\}$ is complete in $L^2(\Omega)$, (4.8) holds for $u = u_j$, and recalling (4.7) we conclude that system (4.1)-(4.2) has a solution. □

5.   Uniform upper bounds

The situation is different from the case without gravity studied in ^[13], where the upper bound could be derived in an elementary way from Hilpert's inequality following ^[15]. Here, we propose a hysteresis variant of the Moser iteration procedure which is new to our knowledge. The price we pay is that we have to restrict the class of admissible Preisach densities $\phi(x, r, v)$ assuming that they have a sufficiently slow decay at infinity, see Hypothesis 3.6.

We first test (4.1) by ${y} = u_i$ and get, using (4.2),

for all $i\in\{1, \dots, n\}$. Let us define the functions

In terms of the sequence $\xi^r_i(x) = {\mathfrak{p}}_r[\lambda, u]_i$ we have

Choosing in (4.5) $z = 0$ and using the fact that the function $\psi$ in (5.2) is increasing, we obtain in both cases $\xi^r_i\ge \xi^r_{i-1}$, $\xi^r_i\le \xi^r_{i-1}$ the inequalities

and (5.1) yields, by Hypothesis 3.4 and the inequality $(v_i - v_{i-1}) v_i \ge (v_i^2 - v_{i-1}^2)/2$,

for $i \in \{1, \dots, n\}$, with a constant $C > 0$ independent of $\tau$. Summing up over $i$ and exploiting the fact that, by the definition of $\xi^r_0$ in (4.6) and the assumptions on $\phi$ and $\lambda$ in Definition 3.1, we have

and we are assuming $v_0 \in L^\infty(\Omega)$, we get the estimate

with a constant $C > 0$ independent of $\tau$.

For parameters $R > 1$ and $k > 1$ to be specified later, we define the function

The function ${U_{k, R}}$ is odd, increasing, continuously differentiable, and its derivative has the form

We test (4.1) with ${y} = {U_{k, R}}(u_i)$, which is an admissible test function, and we obtain

For every $\alpha, \beta \in {\mathbb{R}}$ and for every nondecreasing function $h:{\mathbb{R}} \to {\mathbb{R}}$, we have the elementary inequality

Formula (5.11) for $\alpha = (v_i - v_{i-1})/\tau$, $\beta = v_i$, $h = {U_{k, R}}$ together with (4.2) enables us to estimate

where we put

The next term in (5.10) can be estimate from below by means of Young's inequality

with positive constants $a, b$ independent of $R$ and $k$. To estimate the boundary term in (5.10) from below, we first notice that for all $u\in {\mathbb{R}}$ we have, by virtue of (5.8),

Since $u^*$ is bounded by Hypothesis 3.4, we may use (5.15) and Young's inequality with conjugate exponents $2k+2$ and $(2k+2)/(2k+1)$ to check that there exists a constant $C > 0$ independent of $R$ and $k$ such that

for all $i\in\{1, \dots, n\}$.

Finally, in order to estimate the time-discrete hysteresis term in (5.10), we define the function

Similarly as in (5.4), we have

Summing up in (5.10) over $i = 1, \dots, n$ and using the above estimates (5.12), (5.14), (5.16), and (5.18), we get, using again an estimate of the initial step at $i = 0$ similar to (5.6) and the assumption $v_0 \in L^\infty(\Omega)$,

with constants $a > 0$ and $C > 0$ independent of $\tau$, $k$, and $R$.

For $u\in {\mathbb{R}}$ put

We have

Then

We define an equivalent norm of an element $w$ from the Sobolev space $W^{1, 2}(\Omega)$ by the formula

and using (5.23)-(5.24) rewrite inequality (5.19) in the form

We define

With this notation, by (5.9) we can rewrite the term at the right-hand side of (5.25) as

Let us now compare the $W^{1, 2}$-norm of ${u_i^{(k, R)}}$ with the $W^{1, 2}$-norm of $V_{k, R}(u_i)$ on the left-hand side of (5.25). We have

and, by Young's inequality with conjugate exponents $k$ and $k/(k-1)$,

Hence, recalling (5.21)-(5.22), we estimate

for $i = 1, \dots, n$. From (5.25), (5.27), and (5.7) we thus get

for some constants $C_0, C > 0$ independent of $\tau$, $k$, and $R$. To simplify the notation, we denote by $|\cdot|_q$ the norm in $L^q(\Omega)$. There exists a constant $K > 0$ such that for every $\sigma\in (0, 1)$ the interpolation inequality

holds for every $w \in W^{1, 2}(\Omega)$, see ^[7]. Choosing $\sigma = (C_0 K)^{-1/2}(1+k)^{-1}$, we apply this inequality in (5.29) and obtain

with a constant $C > 0$ depending on $C_0$ and $K$.

So far, this has been a standard Moser argument. The crucial point in the derivation of an upper bound for the time-discrete problem is to find an efficient lower bound for the hysteresis term on the left-hand side of (5.31). By (5.17) we have for $x \in \Omega$ that

where the function ${U_{k, R}}$ is defined in (5.8). Assume first that $u_i(x) > 0$. By (4.5) for all $r > 0$ we have

where $(\cdot)^+$ denotes the positive part. From (5.32) we thus get

In the case $u_i(x) < 0$ we argue similarly and obtain the formula which is valid for both cases

according to Hypothesis 3.6. Our goal is to show that the function

of the argument $u\ge 0$, where we have first used Fubini's theorem and then substituted $z = r+v$, is dominant for $k > (m-3)/2$ over the function

that is, there exists a constant $C_k$ independent of $R$ and possibly dependent on $k$ such that

for all $u \ge 0$. We have in particular

We obviously have (5.36) for $u \le 1$. Indeed, in this case $u \le 1 < R$, and we can compute

so that (5.36) holds with $C_k = 6(k+1)^2$. For $u > 1$ it suffices to prove that

Indeed, then (5.36) follows from the formula $B_{k, R}(u) - B_{k, R}(1) \le C_k (A_{k, R}(u) - A_{k, R}(1))$ and the fact that (5.36) holds for $u = 1$ from the previous step. Note that (5.37) holds automatically for $u > R$. For $u \in (1, R)$ we have by (5.8) that

hence (5.37) and (5.36) hold with $C_k$ as above. Therefore, combining (5.34) and (5.36), in terms of the functions $w_i(x)$ defined in (5.26), estimate (5.31) implies

for all $k > (m-3)/2$ with a constant $C > 0$ independent of $k$ and $R$.

We are ready to start the Moser iterations. Note that we can assume $m\ge3$. Indeed, if (3.16) holds with $m < 3$, it will certainly be true with $m = 3$. Since we have an anisotropic norm on the right-hand side of (5.38), we have to replace the left-hand side by an expression which is compatible with the right-hand side. We shall see that, since

the right choice is

with a constant $C_T$ depending on $T$. We now choose $k_0 > m-3$ and define a sequence $k_j$ for $j \in {\mathbb{N}}$ recurrently by the formula $k_{j} = 2k_{j-1} + 3 - m$, that is,

From (5.39) we obtain

with exponent

Note that $(m-3)/2k_{j-1} \in [0, 1/2)$, so that $\omega_j \in (1/2, 1]$. We define the quantities

Then from (5.41) we get the recurrent relation

for all $j \in {\mathbb{N}}$. We rewrite (5.42) in the form

Summing up the above inequality over $j = 2, \dots, k$ yields

To estimate $L_1$, we observe that by virtue of (5.39) we have

Poincaré's inequality implies

so that from (5.25) and (5.27) it follows that for all $k\ge 1$ we have

By induction over $k$ we prove that for $k_0 > m-3$ as above the estimate

holds true with a constant $C$ depending only on $m$ in Hypothesis 3.6 and on the data of the problem. Hence, by (5.45) and (5.48) we obtain that $L_1 \le C$. Coming back to (5.44), we note also that

and since

we conclude that $\prod_{i = 1}^\infty 1/\omega_i < \infty$. We additionally have

see [23, Theorem 15.5], so that $\prod_{i = 1}^\infty \omega_i > 0$. To summarize, we have

and we conclude from (5.44) that there exists a constant $L^* > 0$ such that $L_k \le L^*$ for all $k\in {\mathbb{N}}$, hence

for a.e. $x \in \Omega$ and all $i\in \{1, \dots, n\}$. By comparison in Eq (4.2), we also have

for a.e. $x \in \Omega$ and all $i\in \{1, \dots, n\}$, with a constant $C > 0$ depending on $U$. Estimate (5.49) additionally implies $|\xi_i^r(x)| \le (\hat{U}-r)^+$ for $\hat{U}{ :=}\max\{U, \Lambda\}$, and we can write in (5.3)

Hence, even if we do not assume the a priori boundedness of $G$ as in (3.9), by assumption (3.7) we obtain

with a constant $C > 0$ independent of $i$ and $\tau$.

6.   Convexity estimate

Recall that the operator $G$ is convexifiable in the sense of Definition 3.3, that is, for every $U > 0$ there exists a twice continuously differentiable mapping $g:[-U, U] \to [-U, U]$ such that $g(0) = 0$, $0 < g_* \le g'(u) \le g^* < \infty$, $|g''(U)| \le \bar{g}$, and $G$ is of the form

where $P$ is a uniformly counterclockwise convex Preisach operator on $[-U, U]$. Let us fix $U$ from (5.49) and the corresponding function $g$. The following result is a straightforward consequence of [12, Proposition 3.6].

Proposition 6.1. Let $P$ be uniformly counterclockwise convex on $[-U, U]$, and let $f$ be an odd increasing function such that $f(0) = 0$. Then there exists $\beta > 0$ such that for every sequence $\{w_i: i = -1, 0, \dots, n-1\}$ in $[-U, U]$ we have

where we set for $w\in {\mathbb{R}}$

We need to define backward steps $u_{-1}$ and $v_{-1}$ satisfying the strong formulation of (4.1)-(4.2) for $i = 0$, that is,

with boundary condition (3.15). Repeating the argument of [12, Proposition 3.3], we use assumptions (3.7) and (3.14) to find for each $0 < \tau < \phi_0(U)/2L^2$ functions $u_{-1}$ and $G[u]_{-1}$ satisfying

as well as, thanks to (3.13), the estimate

with a constant $C > 0$ independent of $\tau$ and $x$. Finally, we construct $v_{-1}$ according to (6.5).

We extend the discrete system (4.1)-(4.2) to $i = 0$ and write it in the form

with $w_i = g(u_i)$, $\theta_i = G[u]_i$, for $i\in \{0, 1, \dots, n\}$ and for an arbitrary test function ${y} \in W^{1, 2}(\Omega)$. We proceed as in ^[12] and test the difference of (6.8) taken at discrete times $i+1$ and $i$

by ${y} = f(w_{i+1} - w_i)$ with

In agreement with (6.3), we have

The hysteresis term is estimated from below by virtue of Proposition 6.1 as follows:

Note that for $|w| \ge \tau({{{\rm{e}}}}^2 - 1)$, we have

so that

We denote by $J$ the set of all $i \in \{0, 1, \dots, n-1\}$ such that $|w_{i+1} - w_i| \ge \tau({{{\rm{e}}}}^2 - 1)$, and by $J^\perp$ its complement. We thus have

hence

with a constant $C > 0$ independent of $x$ and $\tau$. Moreover, for $w_0 \neq w_{-1}$ we have

which yields, together with identity (6.6) and assumption (3.13),

with a constant $C > 0$ independent of $x$ and $\tau$. For $w_0 = w_{-1}$, we interpret $(P[w]_0-P[w]_{-1})/(w_0-w_{-1})$ as $B'_+(w_{-1})$ or $B'_-(w_{-1})$, according to the notation for the Preisach branches introduced after Proposition 3.2, see ^[12] for more details. From (6.14)–(6.17) we thus get

with a constant $C > 0$ independent of $x$ and $\tau$.

We further have by (6.9) and by the monotonicity of $f$ and $g$ that

where in the last step we employed estimates (5.49)-(5.50), and the constant $C > 0$ depends on $U$ from (5.49).

The boundary source term

is bounded by a constant by virtue of Hypothesis 3.4, while the boundary term on the left-hand side

is non-negative by monotonicity of both functions $f$ and $g$. Using (6.18)-(6.19) we thus get from (6.10)

with a constant $C > 0$ independent of $\tau$. We further have

The functions $\kappa$ and $g'$ are bounded and Lipschitz continuous, and

Moreover, since $\theta_i = P[w]_i$ admits a representation similar to (4.3), by Hypothesis 3.4, estimate (5.51), and the Lipschitz continuity of the time-discrete play implied by (4.5), we obtain

whereas, from assumption (3.10),

Thus, employing the Cauchy-Schwarz and Young's inequalities and assumption (3.10) on $g$, we conclude that there exist constants $\delta > 0$ and $K > 0$ independent of $\tau$ such that

The remaining term $(\kappa(x, \theta_{i{+}1}){-} \kappa(x, \theta_{i})){\boldsymbol{\nu}}\cdot \nabla f(w_{i+1} - w_i)$ in (6.20) can be treated as follows. By (3.10) we have

and employing (6.22) and arguing as in (6.24) we get

with some constant $C > 0$ independent of $\tau$. As a consequence of (3.10), (6.20)–(6.25), and (5.7), we thus have the crucial estimate

with a constant $C > 0$ independent of $\tau$.

7.   Proof of Theorem 3.7

The estimates we have derived in Sections 5 and 6 are sufficient for passing to the limit as $\tau\to 0$. For the time-discrete sequence $\{u_i\}$, we repeat literally the compactness argument in Sobolev and Orlicz spaces for piecewise linear and piecewise constant interpolates from [13, Sections 5 and 6]. The most delicate step is the convergence of the hysteresis terms. Estimates (5.7), (5.49), and (6.26) guarantee that the piecewise linear interpolates are bounded in the space

where

and

From [13, Proposition 6.1], the space $\mathcal{B}$ is compactly embedded in $L^1(\Omega; C[0, T])$, and this is enough to ensure the convergence of the hysteresis terms in view of Proposition 3.2. As for the piecewise linear and piecewise constant interpolates of the sequence $\{v_i\}$, their convergence follows from the (strong) convergence of the interpolates of $\{u_i\}$ (compare in (4.2)) combined with the Lebesgue dominated convergence theorem. The limits are the desired solution to Problem (2.24)–(2.26), and this concludes the proof of Theorem 3.7.

8.   Conclusions

In this paper, we have investigated a model for unsaturated flow in a viscoelastic porous medium exhibiting degenerate hysteresis in the pressure-saturation relationship, coupled with gravity-driven moisture flux. This extends previous work on degenerate diffusion in porous media with hysteresis by incorporating the effects of solid matrix deformation and gravity, which introduce significant challenges in the analysis. The key difficulty lies in obtaining a priori bounds for the solutions, which we address by employing a hysteresis variant of the Moser iteration technique. Specifically, we require a condition on the Preisach density, namely, that it decays sufficiently slowly at infinity (Hypothesis 3.6). This condition allows us to establish a uniform upper bound for the solutions. Building upon this estimate, we leverage the convexity and compactness arguments developed in ^[13] to prove the existence of a weak solution to the coupled system. Our approach involves a time discretization scheme, the derivation of estimates independent of the time step, and a passage to the limit.

The results presented here contribute to a deeper understanding of flow phenomena in deformable porous media, relevant to various applications in hydrogeology and related fields. Future work could explore the impact of more complex constitutive models for the solid matrix, as well as the development of efficient numerical methods for the computation of solutions.

Use of Generative-AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The support from the European Union's Horizon Europe research and innovation programme under the Marie Skłodowska-Curie grant agreement No 101102708, from the MŠMT grant 8X23001, and from the GAČR project 24-10586S is gratefully acknowledged.

Conflict of interest

The authors declare no conflicts of interest.