Thierno Sarr^*, Souleye Kane
École Doctorale de Mathématiques et Informatique (ED-MI), Université Cheikh Anta Diop de Dakar, Dakar, Sénégal.
DOI: 10.4236/jamp.2026.146105   PDF    HTML   XML   2 Downloads   45 Views   Citations

Abstract

We study the existence and uniqueness of renormalized solutions for a nonlinear degenerate elliptic-parabolic problem of the form $\frac{?M\left(u\right)}{?t}?\text{div}\left(?u+K\left(u\right)e_z\right)=f\text{?}\text{?}\text{in}\text{?}{Q}_T=\left(0,T\right)\times \Omega ,$ with $L^1$ data $f\in L^1\left(Q_T\right)$ and $u_0\in L^1\left(\Omega \right)$ , arising from the Richards model of unsaturated flow in porous media. Because classical weak solutions fail to be unique at this low-regularity level, we work within the framework of renormalized solutions, which select the physically relevant solution through an energy-dissipation condition at infinity. Existence is established via a double approximation scheme $\left({\psi }_{m,n}\right)$ combined with nonlinear semigroup theory, and uniqueness follows from a comparison principle obtained by the doubling-of-variables technique. We further present a conforming finite-element /implicit-Euler/Picard-iteration scheme validated on the Brooks-Corey soil model, for which we prove an $O\left(h^2+\tau \right)$ error estimate in $L^2\left(\Omega \right)$ and unconditional stability.

Keywords

Renormalized Solutions, Elliptic-Parabolic Equation, Porous Media, Richards Equation, Finite Elements, Picard Iteration, Data, Semigroup Theory

Share and Cite:

1. Introduction

The movement of water in variably saturated soils is governed by the Richards equation [1], which, after the Kirchhoff transformation $u={\displaystyle {\int }_0^h{K}_r\left(\theta \left(s\right)\right)\text{d}s}$ , takes the form

$\frac{\partial M\left(u\right)}{\partial t}-\text{div}\left(\nabla u+K\left(u\right)e_z\right)=f\text{in}{Q}_T=\left(0,T\right)\times \Omega ,$ (1)

supplemented with homogeneous Dirichlet boundary conditions and initial datum $M\left(u\right)\left(\cdot ,0\right)=v_0$ in $\Omega \subset {ℝ}^3$ bounded. Here $M:ℝ\to ℝ$ is continuous, non-decreasing with $M\left(0\right)=0$ , and $K:ℝ\to ℝ$ is Lipschitz with $K\left(0\right)=0$ .

When $f\in L^1\left(Q_T\right)$ and $v_0\in L^1\left(\Omega \right)$ , classical weak solutions are not unique in general [2]. Three equivalent notions have been introduced to restore uniqueness: SOLA [3], entropy solutions [4], and renormalized solutions [5]. The latter, first used for the Boltzmann equation [5] and later adapted to elliptic-parabolic problems by Ammar-Wittbold [6] and Carrillo-Wittbold [7], provide the most flexible framework for handling degenerate diffusions with rough data.

Contributions of this paper.

• We prove existence (Theorem 3.4) and uniqueness (Theorem 3.5) of a renormalized solution of (1) with $L^1$ data, via a two-parameter approximation ${\psi }_{m,n}\left(r\right)=\frac{1}{m}{r}^{+}-\frac{1}{n}{r}^{-}$ and nonlinear semigroup theory (Section 3).

• We design and analyse a finite-element/implicit-Euler/Picard scheme (Section 4), establish an $O\left(h^2+\tau \right)$ error estimate, and validate it on a Brooks-Corey soil model.

Throughout the paper we work under the following standing assumptions:

• $M:ℝ\to ℝ$ is continuous, non-decreasing, $M\left(0\right)=0$ .

• $K:ℝ\to ℝ$ is Lipschitz continuous, $K\left(0\right)=0$ .

• $f\in L^1\left(Q_T\right)$ , $u_0\in L^1\left(\Omega \right)$ .

2. Functional Setting and Preliminary Results

We write $Q_T=\left(0,T\right)\times \Omega$ and denote by $T_k\left(s\right)=\text{max}\left(-k,\text{min}\left(s,k\right)\right)$ the truncation at level $k>0$ , by $H_0$ the Heaviside function, and by sign the sign function. All integrals without explicit domain are over $Q_T$ unless stated otherwise.

2.1. Truncation Operators and Gradient Lemma

Lemma 2.1 (Gradient of truncations, [8]). Let $u$ be measurable on Ω with $T_k\left(u\right)\in W_0^{1,p}\left(\Omega \right)$ for every $k>0$ . Then there exists a unique measurable function $v:\Omega \to {ℝ}^N$ such that

$\nabla T_k\left(u\right)=v{1}_{\left\{\left|u\right|<k\right\}}$ a.e., $\forall k>0.$

Moreover, if $u\in W_0^{1,1}\left(\Omega \right)$ then $v=\nabla u$ .

2.2. Integration-by-Parts Formula

Lemma 2.2 ([9] [10]). Let $M:ℝ\to ℝ$ be continuous and monotone with $M\left(0\right)=0$ . Suppose $u\in L^2\left(0,T;H_0^1\left(\Omega \right)\right)$ with $M\left(u\right)\in L^1\left(Q_T\right)$ , ${\partial }_{t}M\left(u\right)\in L^2\left(0,T;H^{-1}\left(\Omega \right)\right)+L^1\left(Q_T\right)$ , and $u_0\in L^1\left(\Omega \right)$ . Then for every $h\in C_c^{\infty }\left(ℝ\right)$ with $h\left(0\right)=0$ , $\eta \in H_0^1\left(\Omega \right)$ , and $\phi \in C_c^{\infty }\left(Q_T\right)$ with $\phi \left(T\right)=0$ ,

$-{\displaystyle {\int }_0^T\langle \frac{\partial M\left(u\right)}{\partial t},h\left(u-\eta \right)\phi \rangle \text{d}t}={\displaystyle {\int }_{Q_T}\frac{\partial \phi }{\partial t}}{\displaystyle {\int }_u^{u_0}h}\left(r-\eta \right)\text{d}M\left(r\right)\text{d}x\text{d}t.$

2.3. Perturbed Problem and Semigroup Approximation

For $\psi :ℝ\to ℝ$ continuous and strictly increasing with $\psi \left(0\right)=0$ , and $M_k\left(r\right)=M\left(r\right)+\frac{1}{k}r$ (strictly increasing), define the operator

$A_{M_k}^{\psi }:=\left\{\left(M_k\left(u\right),-\text{div}a\left(u,\nabla u\right)+\psi \left(u\right)\right)\right\}$

on the domain $\left\{u\in H_0^1\left(\Omega \right)\cap L^{\infty }\left(\Omega \right):\text{div}a\left(u,\nabla u\right)\in L^{\infty }\left(\Omega \right)\right\}$ , where $a\left(u,\nabla u\right)=\nabla u+K\left(u\right)e_z$ .

Proposition 2.3 ([6] [8]). Under (H1)-(H2), $A_{M_k}^{\psi }$ is $m$ -accretive and densely defined in $L^1\left(\Omega \right)$ . Moreover, as $M_k\to M$ uniformly on compacts, $A_M^{\psi }\subset {\underset{\_}{\lim }}_{k\to \infty }{A}_{M_k}^{\psi }$ .

By nonlinear semigroup theory [11], Proposition 2.3 yields a unique mild solution $v_k=M_k\left(u_k\right)\in C\left(\left[0,T\right];L^1\left(\Omega \right)\right)$ of the Cauchy problem ${\dot{v}}_k+A_{M_k}^{\psi }\left(v_k\right)\ni f$ , $v_k\left(0\right)=v_{0k}$ . As $k\to \infty$ and $v_{0k}\to v_0$ in $L^1\left(\Omega \right)$ , we have $v_k\to v=M\left(u\right)$ in $C\left(\left[0,T\right];L^1\left(\Omega \right)\right)$ .

3. Renormalized Solutions: Existence and Uniqueness

3.1. Definition

Definition 3.1 (Renormalized solution). A measurable function $u:Q_T\to ℝ$ is a renormalized solution of (1) if

(i) $M\left(u\right)\in L^1\left(Q_T\right)$ ,

(ii) $T_k\left(u\right)\in L^2\left(0,T;H_0^1\left(\Omega \right)\right)$ for every $k>0$ ,

(iii) for every $h\in C_c^1\left(ℝ\right)$ and $\phi \in C_c^{\infty }\left(\left(0,T\right)\times \Omega \right)$ ,

$\begin{array}{l}{\displaystyle {\int }_{Q_T}\frac{\partial \phi }{\partial t}}{\displaystyle {\int }_u^{u_0}h}\left(r\right)\text{d}M\left(r\right)\text{d}x\text{d}t+{\displaystyle {\int }_{Q_T}fh\left(u\right)\phi \text{d}x\text{d}t}\\ ={\displaystyle {\int }_{Q_T}\left(\nabla u+K\left(u\right)e_z\right)\cdot \nabla \left(h\left(u\right)\phi \right)\text{d}x\text{d}t},\end{array}$ (2)

(iv) the following energy condition at infinity holds:

${\displaystyle {\int }_{Q_T\cap \left\{k\le \left|u\right|\le k+1\right\}}\left(\nabla u+K\left(u\right)e_z\right)\cdot \nabla u\text{d}x\text{d}t}\stackrel{k\to +\infty }{\to }0.$(3)

Remark 3.2. Condition (3) selects the physical solution among (possibly infinitely many) weak solutions: it expresses that no energy is dissipated at infinite amplitude levels. When $\text{supp}h\subset \left[-k,k\right]$ , the term $\left(\nabla u+K\left(u\right)e_z\right)\cdot \nabla \left(h\left(u\right)\phi \right)$ is identified with $\left(\nabla T_k\left(u\right)+K\left(T_k\left(u\right)\right)e_z\right)\cdot \nabla \left(h\left(T_k\left(u\right)\right)\phi \right)$ , which is well-defined by (ii).

3.2. Comparison Lemma

Lemma 3.3 (Monotone comparison). Let $u_0,{\tilde{u}}_0,f,\tilde{f}\in L^{\infty }\left(Q_T\right)$ and let $\psi ,\tilde{\psi }:ℝ\to ℝ$ be continuous and strictly increasing with $\psi \left(0\right)=\tilde{\psi }\left(0\right)=0$ . Let $u,\tilde{u}$ be weak solutions of the perturbed problem (1) with perturbations $\psi$ and $\tilde{\psi }$ respectively. If $u_0\le {\tilde{u}}_0$ a.e. on Ω, $f\le \tilde{f}$ a.e. on $Q_T$ , and $\tilde{\psi }\left(r\right)\le \psi \left(r\right)$ for all $r\in ℝ$ , then $u\le \tilde{u}$ a.e. on $Q_T$ .

Proof. Taking $\frac{1}{L}\phi \left(t\right){\rho }_p\left(t-s\right)T_L\left(u\left(t,x\right)-\tilde{u}\left(s,x\right)\right)$ as test function and using the identity

$\left(\psi \left(u\right)-\tilde{\psi }\left(\tilde{u}\right)\right)\text{sign}\left(u-\tilde{u}\right)\ge {\left(\tilde{\psi }\left(u\right)-\tilde{\psi }\left(\tilde{u}\right)\right)}^{+},$

together with the same doubling argument as in Lemma 2.2 and Proposition 2.3, one obtains ${\displaystyle {\int }_{Q_T}{\left(\tilde{\psi }\left(u\right)-\tilde{\psi }\left(\tilde{u}\right)\right)}^{+}\text{d}x\text{d}t}\le 0$. Since $\tilde{\psi }$ is strictly increasing, $u\le \tilde{u}$ a.e. □

3.3. Existence

Theorem 3.4 (Existence). Under (H1)-(H3), there exists at least one renormalized solution of (1).

Proof. Step 1: Approximation. Set $f_{m,n}=\text{max}\left(\text{min}\left(f,m\right),-n\right)$ , $u_{0,m,n}=\text{max}\left(\text{min}\left(u_0,m\right),-n\right)$ , and ${\psi }_{m,n}\left(r\right)=\frac{1}{m}{r}^{+}-\frac{1}{n}{r}^{-}$ for $m,n\in {\mathbb{N}}^{\text{*}}$ . By Proposition 2.3 and nonlinear semigroup theory (applied with $\psi ={\psi }_{m,n}$ , which is strictly increasing), there exists a unique weak solution $u_{m,n}$ of the perturbed problem

$\frac{\partial M\left(u\right)}{\partial t}+{\psi }_{m,n}\left(u\right)-\text{div}a\left(u,\nabla u\right)=f_{m,n}\text{in}{Q}_T,M\left(u\right)\left(\cdot ,0\right)=u_{0,m,n},{u|}_{\partial \Omega }=0.$ (4)

Step 2: Monotonicity and pointwise limits. By Lemma 3.3,

$u_{m,n}{\uparrow }_m{u}_n\text{and}{u}_n{\downarrow }_{n}u\text{a}.\text{e}.\text{on}{Q}_T,$

where $u_n$ and $u$ are measurable and almost-everywhere finite (see [6]).

Step 3: Energy estimate. Choosing test function $T_k\left(u_{m,n}\right){\phi }_{\delta }$ in (4), where ${\phi }_{\delta }\in C^1\left(\left[0,T\right]\right)$ is a cut-off near $T$ , using Lemma 2.2 and the monotonicity of ${\psi }_{m,n}$ , one obtains

$C{\displaystyle {\int }_0^T{\displaystyle {\int }_{\Omega }{\left|\nabla T_k\left(u_{m,n}\right)\right|}^2\text{d}x\text{d}t}}\le k\left({\Vert f\Vert }_{L^1\left(Q_T\right)}+{\Vert u_0\Vert }_{L^1\left(\Omega \right)}\right).$ (5)

By the Poincaré inequality, $\left(T_k\left(u_{m,n}\right)\right)$ is bounded in $L^2\left(0,T;H_0^1\left(\Omega \right)\right)$ , so $T_k\left(u_{m,n}\right)\rightharpoonup T_k\left(u\right)$ weakly in $L^2\left(0,T;H_0^1\left(\Omega \right)\right)$ .

Step 4: Strong convergence of truncations. Taking $h\left(u_n\right)\left(T_k\left(u_n\right)-T_k\left(u\right)\right)\phi$ as test function in the equation satisfied by $u_n$ , one shows that $I_1+I_2+I_3\to 0$ (where $I_1$ involves the time derivative, $I_2$ the gradient cross-terms, and $I_3$ the convection), yielding ${\Vert T_k\left(u_n\right)\Vert }_{L^2\left(H_0^1\right)}\to {\Vert T_k\left(u\right)\Vert }_{L^2\left(H_0^1\right)}$ . Combined with the weak convergence, this gives

$T_k\left(u_n\right)\to T_k\left(u\right)\text{strongly}\text{in}{L}^2\left(0,T;H_0^1\left(\Omega \right)\right).$ (6)

Step 5: Passage to the limit. Using the strong convergence (6) and the Lebesgue dominated convergence theorem, one passes to the limit in each term of the equation for $u_n$ as $n\to \infty$ to obtain (2).

Step 6: Energy condition at infinity. Taking $\left(T_{k+1}\left(u_n\right)-T_k\left(u_n\right)\right)\phi$ as test function and letting $n\to \infty$ , then $k\to \infty$ with $\phi =1$ , yields (3). □

3.4. Uniqueness

Theorem 3.5 (Uniqueness and comparison). Under (H1)-(H3), for $i=1,2$ let $u_{0i}\in L^1\left(\Omega \right)$ , $f_i\in L^1\left(Q_T\right)$ , and $u_i$ be renormalized solutions of (1) with data $\left(u_{0i},f_i\right)$ . Then, for all $t\in \left(0,T\right)$ ,

${\displaystyle {\int }_{\Omega }{\left(M\left(u_1\right)\left(t\right)-M\left(u_2\right)\left(t\right)\right)}^{+}\text{d}x}\le {\displaystyle {\int }_{\Omega }{\left(u_{01}-u_{02}\right)}^{+}\text{d}x}+{\displaystyle {\int }_0^t{\displaystyle {\int }_{\Omega }{\left(f_1-f_2\right)}^{+}\text{d}x\text{d}t}}.$ (7)

In particular, $u_1=u_2$ a.e. when $u_{01}=u_{02}$ and $f_1=f_2$ .

Proof sketch. The proof uses the double-variable technique of [7], exploiting the structural condition

$\left(a\left(r,\zeta \right)-a\left(s,\eta \right)\right)\cdot \left(\zeta -\eta \right)+C\left(r,s\right)\left(1+{\left|\zeta \right|}^p+{\left|\eta \right|}^p\right)\left|r-s\right|\ge \Gamma \left(r,s\right)\cdot \zeta +\widehat{\Gamma }\left(r,s\right)\cdot \eta ,$

which holds for $a\left(u,\nabla u\right)=\nabla u+K\left(u\right)e_z$ since $K$ is Lipschitz. A doubling of both the space and time variables, combined with the energy condition (3), yields (7). □

4. Numerical Scheme and Analysis

Since the domain is one-dimensional, $\Omega =\left]0,L\right[\subset ℝ$ , the vertical coordinate $z$ of the analytical model coincides with the spatial coordinate $x$ used in the finite-element discretisation. Consequently, the gravitational convection term $K\left(u\right){\partial }_z$ appearing in the flux $\nabla u+K\left(u\right)e_z$ is written in the finite-element coordinates as

$K\left(u\right){\partial }_z\equiv K\left(u\right)\frac{\partial }{\partial x},$

so that the weak contribution of this term against a test function ${\phi }_l$ reads ${\displaystyle {\int }_{\Omega }K\left(u_h\right)\frac{\partial {\phi }_l}{\partial x}\text{d}x}$ .

4.1. Spatial Discretisation by Finite Elements

Let $\Omega =\left]0,L\right[$ with $L=1\text{m}$ . We construct a mesh ${\mathcal{T}}_h=\left\{L_1,\cdots ,L_{N+1}\right\}$ with nodes $0=x_0<x_1<\cdots <x_{N+1}=L$ and $h={\text{max}}_i\left(x_{i+1}-x_i\right)$ . The conforming finite-element space is

$V_h=\left\{v_h\in C^0\left(\overline{\Omega }\right):{v_h|}_{L_i}\in 1_{},v_h\left(0\right)=v_h\left(L\right)=0\right\},$

with nodal basis ${\left\{{\phi }_j\right\}}_{j=1}^N$ . The semi-discrete problem reads: find $u_h\left(t\right)={\displaystyle {\sum }_{m=1}^N{u}_m\left(t\right){\phi }_m}\in V_h$ such that

${\displaystyle {\int }_{\Omega }\frac{\partial M\left(u_h\right)}{\partial t}{\phi }_l\text{d}x}+{\displaystyle {\int }_{\Omega }\nabla u_h\cdot \nabla {\phi }_l\text{d}x}+{\displaystyle {\int }_{\Omega }K\left(u_h\right){\partial }_z{\phi }_l\text{d}x}={\displaystyle {\int }_{\Omega }f{\phi }_l\text{d}x},l=1,\cdots ,N.$ (8)

4.2. Temporal Discretisation and Picard Linearisation

With time step $\tau =T/N_t$ and $u_h^k\approx u_h\left(\cdot ,t_k\right)$ , the implicit Euler step gives:

${\displaystyle {\int }_{\Omega }\frac{M\left(u^k\right)-M\left(u^{k-1}\right)}{\tau }{\phi }_l\text{d}x}+{\displaystyle {\int }_{\Omega }\nabla u^k\cdot \nabla {\phi }_l\text{d}x}+{\displaystyle {\int }_{\Omega }K\left(u^k\right){\partial }_z{\phi }_l\text{d}x}={\displaystyle {\int }_{\Omega }{f}^k{\phi }_l\text{d}x}.$ (9)

Since (9) is nonlinear in $u^k$ , we linearise via a first-order Taylor expansion (Picard method):

$M\left(u^{k,i}\right)\approx M\left(u^{k,i-1}\right)+C_1^{k,i-1}\left(u^{k,i}-u^{k,i-1}\right),K\left(u^{k,i}\right)\approx K\left(u^{k,i-1}\right)+C_2^{k,i-1}\left(u^{k,i}-u^{k,i-1}\right),$

where $C_j^{k,i-1}={{\partial }_u\left\{M,K\right\}|}_{u^{k,i-1}}$ . This yields, at each Picard step $i$ , the linear system $A^{k,i-1}{u}^{k,i}=B^{k,i-1}$ with stiffness matrix

$A_{lm}^{k,i-1}={\displaystyle {\int }_{\Omega }{C}_1^{k,i-1}{\phi }_l{\phi }_m\text{d}x}+\tau {\displaystyle {\int }_{\Omega }\nabla {\phi }_l\cdot \nabla {\phi }_m\text{d}x}+\tau {\displaystyle {\int }_{\Omega }{C}_2^{k,i-1}{\phi }_l{\partial }_z{\phi }_m\text{d}x}.$

4.3. Error Analysis

We work under the additional regularity hypotheses:

• $M$ is Lipschitz with $M^{\prime }\ge m_0>0$ on every compact (non-degenerate regime of Brooks-Corey [12]).

• $u\in L^{\infty }\left(0,T;H^2\left(\Omega \right)\right)\cap H^1\left(0,T;H_0^1\left(\Omega \right)\right)\cap H^2\left(0,T;L^2\left(\Omega \right)\right)$ .

• The Picard iteration converges in a uniformly bounded number $N_{\text{Picard}}$ of steps (guaranteed when $\tau$ is small relative to the Lipschitz constant of $M^{\prime }$ ).

Theorem 4.1 (Convergence order). Under (H1)-(H6), there exists $C>0$ independent of $h$ and $\tau$ such that

$\underset{0\le k\le N_t}{\max }{\Vert u\left(\cdot ,t_k\right)-u_h^k\Vert }_{L^2\left(\Omega \right)}\le C\left(h^2+\tau \right),$ (10)

${\left(\tau {\displaystyle \sum }_{k=1}^{N_t}{\Vert u\left(\cdot ,t_k\right)-u_h^k\Vert }_{H_0^1\left(\text{Ω}\right)}^2\right)}^{1/2}\le C\left(h+\tau \right).$ (11)

Proof sketch. Decompose $u\left(\cdot ,t_k\right)-u_h^k={\eta }^k+{\theta }^k$ where ${\eta }^k=u\left(\cdot ,t_k\right)-{\Pi }_{h}u\left(\cdot ,t_k\right)$ is the elliptic projection error and ${\theta }^k={\Pi }_{h}u\left(\cdot ,t_k\right)-u_h^k$ . By standard interpolation theory and $H^2$ -regularity of $\Omega =\left]0,L\right[$ ,

${\Vert {\eta }^k\Vert }_{L^2}\le C{h}^2{\Vert u\left(\cdot ,t_k\right)\Vert }_{H^2},{\Vert {\eta }^k\Vert }_{H^1}\le Ch{\Vert u\left(\cdot ,t_k\right)\Vert }_{H^2}.$

Testing the error equation for ${\theta }^k$ with ${\theta }^k$ and using the monotonicity of $M$ (hypothesis (H4)) gives

${\Vert {\theta }^k\Vert }_{L^2}^2\le {\text{e}}^{C{t}_k}{\Vert {\theta }^0\Vert }_{L^2}^2+C{\displaystyle \sum }_{j=1}^k\tau \left({\tau }^2+h^4\right).$

Since ${\theta }^0=0$ , we obtain (10). The estimate (11) follows analogously without the Aubin-Nitsche duality argument. □

Proposition 4.2 (Unconditional stability). Under (H1)-(H2), for every $\tau >0$ and $h>0$ ,

$\underset{0\le k\le N_t}{\max }{\Vert u_h^k\Vert }_{L^2\left(\Omega \right)}^2+\tau {\displaystyle \sum }_{k=1}^{N_t}{\Vert \nabla u_h^k\Vert }_{L^2\left(\Omega \right)}^2\le C\left({\Vert u_0\Vert }_{L^2\left(\Omega \right)}^2+{\Vert f\Vert }_{L^2\left(Q_T\right)}^2\right).$

No CFL condition is required.

5. Numerical Experiment: Brooks-Corey Soil Model

5.1. Physical Parameters

We validate the scheme on a one-dimensional infiltration problem with Brooks-Corey constitutive relations [12] (see Table 1 for the physical parameters and numerical settings):

$M\left(u\right)=\theta \left(h\right)=\{\begin{array}{ll}{\theta }_r+\left({\theta }_s-{\theta }_r\right){\left(\frac{h}{h_c}\right)}^{-\beta }\hfill & \text{if}h<h_c\hfill \\ {\theta }_s\hfill & \text{if}h\ge h_c\hfill \end{array},K_r\left(\theta \right)={\left(\frac{\theta -{\theta }_r}{{\theta }_s-{\theta }_r}\right)}^{3+2/\beta }.$

Table 1. Brooks-Corey parameters and numerical settings.

Parameter	Symbol	Value	Unit
Residual water content	${\theta }_r$	0.064	–
Saturated water content	${\theta }_s$	0.126	–
Air-entry potential	$h_c$	−0.078	m
Pore-size index	$\beta$	24.774	–
Saturated conductivity	$K_{\text{sat}}$	0.0095	m/s
Time step	$\tau$	0.01	s
Mesh size	$h$	0.01	m

5.2. Initial and Boundary Conditions

Reduction to homogeneous Dirichlet data. The analytical model (1) is set with homogeneous Dirichlet boundary conditions ${u|}_{\partial \Omega }=0$ and the finite-element space $V_h$ enforces $v_h\left(0\right)=v_h\left(L\right)=0$ . In the numerical experiment the physical boundary data are $u\left(0,t\right)=0$ and $u\left(1,t\right)=h_c$ . To reconcile these with the analytical framework, we introduce the affine lifting

$\tilde{u}\left(x,t\right)=u\left(x,t\right)-h_{c}x,x\in \left[0,1\right],$

which satisfies $\tilde{u}\left(0,t\right)=0$ and $\tilde{u}\left(1,t\right)=0$ , so that $\tilde{u}$ belongs to the homogeneous space $V_h$ for all $t$ . The equation for $\tilde{u}$ retains the same form as (1) with a modified right-hand side $\tilde{f}=f+\text{Δ}\left(h_{c}x\right)=f$ (since the lifting is affine), thereby preserving the structure of the analytical model. In the remainder of this section, $u$ denotes the lifted variable $\tilde{u}$ and all results of Sections 3-4 apply without modification.

Initial condition. The initial condition is chosen as

$u_0\left(x\right)=h_c\left(1-{\cos }^2\left(\frac{\pi x}{2}\right)\right),x\in \left[0,1\right],$

which satisfies $u_0\left(0\right)=0$ (surface at saturation) and $u_0\left(1\right)=h_c$ (unsaturated base), consistently with the boundary data above. Its derivative ${u^{\prime }}_0\left(x\right)=\frac{\pi }{2}{h}_c\sin \left(\pi x\right)\le 0$ reflects a physically realistic suction profile increasing with depth.

5.3. Results and Discussion

Figure 1 shows the evolution of the matric potential $u\left(x,t\right)$ over ten time steps.

Figure 1. Evolution of the matric potential $u\left(x,t\right)$ over ten time steps ($\tau =0.01\text{s}$ ). The blue curve is the initial condition; as $t\to \infty$ the profile tends uniformly to zero (fully saturated medium). Computed with Scilab on a mesh of $N=100$ nodes.

The simulation exhibits three physically consistent features:

1) The initial profile (blue) displays a strong gradient near $x=1$ , consistent with the prescribed base condition.

2) As $t$ increases, $u$ converges uniformly toward zero, reflecting the progressive saturation of the medium.

3) Convergence is rapid in the early steps and slows as the profile becomes quasi-uniform—a behaviour well-known from Richards equation theory.

Unconditional stability (Proposition 4.2) justifies the choice $\tau =h=0.01$ independently of any CFL constraint.

6. Conclusions

We have established the well-posedness of the Richards-type elliptic-parabolic problem (1) in the renormalized sense for $L^1$ data, and designed a provably convergent and unconditionally stable finite-element scheme. The two-parameter approximation ${\psi }_{m,n}$ provides a flexible tool for handling the degeneracy of $M$ and the lack of integrability of the data simultaneously.

Several extensions are natural:

• Multidimensional simulations on unstructured meshes (FEM or FVM).

• Extension to the fully degenerate case $m_0=0$ (Richards equation at full saturation) using time-space estimates of [10].

• Adaptive time-stepping guided by the Picard convergence rate.

Acknowledgements

The authors thank the members of the Numerical Analysis and PDE Research Group at UCAD for stimulating discussions.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Adams, R.A. and Fournier, J.J.F. (2003) Sobolev Spaces. 2nd Edition, Academic Press.
[2]	Bénilan, P., Boccardo, L., Gallouët, T., Gariepy, R., Pierre, M. and Vázquez, J.L. (1995) An L1-Theory of Existence and Uniqueness of Solutions of Nonlinear Elliptic Equations. Annali della Scuola Normale Superiore di Pisa, 22, 241-273.
[3]	Dall’Aglio, A. (1996) Approximated Solutions of Equations with L1 Data. Application to the H-Convergence of Quasi-Linear Parabolic Equations. Annali di Matematica Pura ed Applicata, 170, 207-240.[CrossRef]
[4]	Carrillo, J. (1999) Entropy Solutions for Nonlinear Degenerate Problems. Archive for Rational Mechanics and Analysis, 147, 269-361.[CrossRef]
[5]	DiPerna, R.J. and Lions, P.L. (1989) On the Cauchy Problem for Boltzmann Equations: Global Existence and Weak Stability. The Annals of Mathematics, 130, 321-366.[CrossRef]
[6]	Ammar, K. and Wittbold, P. (2003) Existence of Renormalized Solutions of Degenerate Elliptic-Parabolic Problems. Proceedings of the Royal Society of Edinburgh: Section A Mathematics, 133, 477-496.[CrossRef]
[7]	Carrillo, J. and Wittbold, P. (1999) Uniqueness of Renormalized Solutions of Degenerate Elliptic-Parabolic Problems. Journal of Differential Equations, 156, 93-121.[CrossRef]
[8]	Benilan, P. and Wittbold, P. (1996) On Mild and Weak Solutions of Elliptic-Parabolic Problems. Advances in Differential Equations, 1, 1053-1073.[CrossRef]
[9]	Wilhelm Alt, H. and Luckhaus, S. (1983) Quasilinear Elliptic-Parabolic Differential Equations. Mathematische Zeitschrift, 183, 311-341.[CrossRef]
[10]	Otto, F. (1996) L1-Contraction and Uniqueness for Quasilinear Elliptic-Parabolic Equations. Journal of Differential Equations, 131, 20-38.[CrossRef]
[11]	Brezis, H. (1983) Analyse fonctionnelle: Théorie et applications. Masson.
[12]	Brooks, R.H. and Corey, A.T. (1966) Properties of Porous Media Affecting Fluid Flow. Journal of the Irrigation and Drainage Division, 92, 61-88.[CrossRef]